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1. INTRODUCTION 


Project Objectives 

The primary objectives of this research project, as stated in the initial proposal by 
the first author of this report, were to develop the following computational tools for use 
in active tip-clearance control: 

• A discrete-layer plate theory (both semi-analytic and finite element approximations) 
for laminated composites structures with embedded piezoelectric layers. 

• A 3-D general discrete- layer element generated in curvilinear coordinates with con- 
stant out-of-plane displacement components for modeling laminated composite piezo- 
electric shells with curved surfaces. This tool will be used to represent the general 
geometries of the engine blades and the casing. 

This report describes the theoretical developments and implementation of these models as 
well as additional work completed as part of this study. Conclusions regarding the nature 
and performance of these models are presented, and additional testing and applications 
are suggested. 

Structure of Report 

The primary components of this research were the theoretical development and intro- 
duction of finite element models and semi-analytic solutions using discrete-layer approxi- 
mations for the analysis of laminated piezoelectric plates and shells. These models form 
the basic thrust and results of the present study. However, as part of this development, 
it was found necessary to develop an additional component as both a check and a com- 
putational alternative for the plate geometry. This additional component involved the 
development of exact solutions for the static and dynamic behavior of simply-supported 
laminated piezoelectric plates. Although not included in the original proposal, the de- 
velopment of the exact solutions was critical to the successful development of the plate 
element because there are effectively no other results in the literature for laminated piezo- 
electric plates other than very simple approximate models. The exact solutions filled a 
needed gap in the understanding of these laminates, and because their development was 
completed under the auspices of this project, the results are included here. 

The following chapter gives a brief review of the literature to describe the current state- 
of-the-art in the field of control using active materials, with emphasis given to piezoelectric 
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materials. Also given in this chapter is an overview of the linear theory of piezoelectricity 
and the governing equations to be used in the following chapters. In the third chap- 
ter, semi-analytic solutions developed for the plate geometry are described. The fourth 
and fifth chapters give the development of the plate and shell elements, respectively, for 
laminated piezoelectric solids using the discrete-layer theory. The sixth chapter details 
the development of the exact solutions for the static and dynamic behavior of simply- 
supported piezoelectric laminates. The seventh chapter contains representative example 
problems for the models developed in the previous chapters. A discussion of potential 
future work and summary are given in the eighth and final chapter. 
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2. BACKGROUND 


This section includes a review of the literature regarding both discrete-layer theories 
and so-called smart structures as these topics apply to the present study. The fundamental 
equations related to the linear theory of piezoelectricity are also included in this section, 
and are provided as a useful foundation for the chapters that follow. 

Literature Review 

. Piezoelectricity 

In 1880, Pierre and Jacques Curie discovered that the electrical charges on certain 
dielectric crystal could be produced by impressing mechanical stress. Later in 1881, 
Hankie named this direct and its converse effect “piezoelectricity”. With the knowledge 
of solid mechanics and electricity, the mathematical foundations and governing equations 
were available and used in the original work of Voight [63]. In addition to this volume, two 
other classical works of Mason [64] and Cady [66], are cited frequently in latter studies of 
piezoelectricity. Both works gave a special attention in the physical properties of crystals, 
as well as practical applications. 

Several works published in 1960’s initiated the numerical solution of problems in- 
volving piezoelectric solids. In 1967, Eer Nisse introduced the variational calculus and 
a Ritz-approximation for analyzing electroelastic vibration problems [5, 6]. In his work, 
variational formulas were derived for the equations of motion of piezoelectric cylinders 
and were applied to find the resonant frequencies of thick discs. This application gave 
a beginning for latter studies using numerical approximation. The monograph of Tier- 
sten [67], provides a comprehensive treatment and application of the linear theory of 
piezoelectricity. In this book, a systematic derivation of the governing equations of linear 
piezoelectricity were given, and results of homogeneous plate vibration problems were 
shown. In later works dealing with piezoelectric vibrations, this monograph is without 
question the most oft-cited reference. 

On the basis of the works of Eer Nisse [5, 6], Holland [7] presented an application to 
rectilinear geometry. With specified trial functions, the resonant properties of piezoelec- 
tric ceramic rectangular parallelepipeds were discussed, as well as several mode shapes 
were contoured. To avoid the limitation of using Ritz-approximation under complicated 
boundary conditions, Allik and Hughes [9] introduced a finite element, formulation for the 
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equations of piezoelectricity. A electroelastic tetrahedral element was also developed in 
this work but without further application and numerical results. This work began the ef- 
forts of using modern computational techniques to represent the behavior of piezoelectric 
solids. 

Several studies have appeared using finite element approximations to the equations 
of linear piezoelectricity. Kagawa and Yamabuchi [14] developed a computer program to 
solve the axisymmetric vibrations of a piezoelectric circular rod of finite length. Later in 
the 1980’s, several piezoelectric structures were analyzed by Naillon et al. [20], as well 
as three-dimensional structures by Ostergaard and Pawlak [25]. Being different from Eer 
Nisse’s approximation [6], Kunkel et al. [38] used the finite element method to find the 
natural vibrational mode of axially symmetric piezoelectric ceramic disks. 

Three volumes, which are not referred here but have been cited frequently in recent 
studies, are the works of Nye [6S], Bottom [71] and Zelenka [?2]. These books provide 
details and more recent information on the properties of (quartz) crystals. In Zelenka’s 
book, Piezoelectric resonators and non-linear properties of crystal are an area of focus. 

Theories of Thick Laminates 

The classical laminated plate theory (CLPT), which is a direct extension of classical 
plate theory (CPT), has been used to solve many composite plate problems for years. 
Because the assumed Kirchhoff hypothesis remains in CLPT, the thickness effect, which 
is known to cause poor results when analyzing thick plates with the use of CPT, similarly 
occurs when using CLPT to analyze thick laminates. According to the literature reviews 
of several studies [32, 43, 44], the initial work to take into account the thickness effect 
could be credited to Basset [1] in 1890. After this work, several displacement-based 
refined theories developed by Hidebrand, Reissner and Thomas [2], as well as Mindlin 
[3], were extremely useful to later studies dealing with the transverse stresses in plates. 
Another commonly used book is that of Timoshenko [65], which gives development of 
a generalized plate theory and some elasticity solutions for specific plate geometries. In 
many later studies to refine the CLPT, the above works are frequently cited. 

In Pagano’s works [10, 11], limitations of the CLPT for laminated elastic plates were 
investigated. Solutions of several specific composite problems using CLPT were com- 
pared with the corresponding elasticity solutions (i.e. exact solutions). The conclusions 
pointed out the undesirable effects and limitations of CLPT, and motivated later re- 
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finment. To deal with the shear deformation and rotatory inertia effects in arbitrarily 
laminated anisotropic plates, the first-order shear deformation plate theory (FSDT) [2, 3] 
was generalized by Yang et al. [4]. In this generalized FSDT, the normals to the mid- 
plane before deformation remains straight (i.e. first-order formulation) but not necessarily 
normal to the mid-plane after deformation. However, a correction to transverse stiffness 
was required to satisfy the boundary conditions (i.e. the zero transverse shear stress on 
the top and bottom surfaces of plate). Several applications using this theory in bending, 
vibration and transient analyses were shown by Whitney et al. [12] and Reddy [17, 19]. 
Following the generalized FSDT, higher order formulations were given to obtain a better 
description for thick laminates. A variationally consistent high order shear deformation 
laminate plate theory was introduced by Reddy [22, 27]. Using this theory, the correction 
needed in the FSDT was not required. Later, for reasons of efficiency, Putcha and Reddy 
developed a mixed shear flexible finite element [27] to apply the high order theories. An 
analysis of stability and natural vibration for laminated plates were also represented. In 
addition to Reddy’s studies, Khdeir presented several applications using his refined shear 
deformation theory [40]. Results available in the literature have been adapted frequently 
for comparisons in many later studies. 

Discrete-layer Theories 

Several works employed to develop so-called discrete-layer theories are specially intro- 
duced here for their usefulness in this research. The most initial work that use a layerwise 
displacement theory to analyze layered piezoelectric plates and layered anisotropic plates 
could be credited to Pauley [13] in 1974. In his dissertation, a layer-wise analysis was 
used to study the free vibration characteristics of infinite laminated piezoelectric plates 
using finite element thickness approximations. However, the analysis was limited under a 
condition of plane strain. Even though this work was a precursor to later development of 
discrete-layer theories and research in laminated piezoelectric structures, no later studies 
have cited this effort. 

In the work of Reddy [32], a general two-dimensional shear deformation theory of 
laminate plates (GLPT) was represented. Based on GLPT, a desired degree of approxi- 
mation of the displacements through the laminate thickness can be given. Additionally, 
this work indicated that the GLPT could represent a generalized form of many other 
laminate theories. The CLPT, the generalized FSDT, the higher order theory [22], and 
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several other refined laminate theories can be considered as special cases deduced from 
GLPT. The discrete-layer theories are also extended form the GLPT. 

After GLPT was developed, Reddy et al. [43, 44] presented several applications and 
detailed discussion. In [43], a plate bending element was developed using GLPT. The 
accuracy of using this element was also evaluated by comparison with the exact solutions 
of the generalized plate theory [65] and 3D-elasticity solutions. When using the bending 
element, a layer-wise description of the inplane displacements and stresses of laminate 
can be given, and an improved approximation to the transverse shear deformation was 
obtained. Results of several examples showed that the thickness effect was effectively 
eliminated. Almost at the same time, another work [44] represented the exact analytic 
solutions of GLPT in two cases, cylindrical bending and simple support plate. By com- 
paring with the 3D-elasticity solutions, this work confirmed that GLPT allowed accurate 
determination of interlaminar stresses. 

In the work of Robbins and Reddy [53], the GLPT was adapted to laminated piezo- 
electric beams. This could be the first work that involved the GLPT into the analysis of 
the smart material. Four different displacement-based finite element models, which were 
derived from GLPT, were represented with numerical results. However, the piezoelectric 
effect was only modeled using an induced strain rather than the equation of piezoelectric- 
ity. 

Finally, the recent work of Heyliger and Saravanos [59, 62] represented discrete-layer 
theories that can deal with the thick piezoelectric laminates. The kinematic assumptions 
were based on the GLPT. To describe the in-plane performance of plates or beams, piece- 
wise linear variations of the components were used through the composite thickness. In 
addition, finite element and global/Ritz approximation in plane of laminate were involved 
to solve several test problems. 

Intelligent Structures And Piezoelectric Laminate Theories 

In this section, the primary characteristics of piezoelectric beam, plate, and shell 
theories are introduced. 

Compared with the use of piezoelectrics in transducer applications, the study of dis- 
tributed piezoelectric actuators for all solids is quite recent, with most papers appearing 
after 1985. The work of Bailey and Hubbard [23] may be the first that introduced dis- 
tributed piezoelectric polymers as actuators to control the bending vibration of cantilever 
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beams. Fanson and Chen [28, 29] also demonstrated that the use of piezoelectric materials 
as sensors/actuators in beam vibration control was feasible. The passive structure ele- 
ment for the control of large space structures can also be replaced by piezoelectric active 
members. 

Later, Crawley and Luis [30] presented an analytic and experimental development of 
one-dimensional piezoelectric elements. These elements can be placed either on the sur- 
face or embedded within structural laminated beams, and function as actuators to excite 
the steady-state resonant vibrations in the cantilevered beams. The given conclusion em- 
phasized that the existence of the embedded actuators may not affect the elastic modular 
of the composite structures, but would reduce the ultimate strength of the laminate by 
20%. Finally, the work of Robbins and Reddy [53] is repeated here for their analysis of 
simulated piezoelectric laminated beams. 

Beginning in 1987, a series of publications by Lee et al. [39, 84. 47] initiated efforts 
to develop theoretical models for bending and torsional control in laminated piezoelectric 
plates. In [39], Lee and Moon presented a set of piezopolymer devices for the model control 
of piezoelectric laminates. Experimental results were followed and used to compare with 
the theoretical predictions. The latter work by Lee [84] introduced in detail the governing 
equations for piezoelectric laminate plates, as well as the reciprocal relationships of the 
piezoelectric sensor/actuator. The assumptions of CLPT were involved because only 
slender plates were considered. 

Wang and Rogers [52] also used the assumptions of CLPT to model laminated plates 
with spatially distributed piezoelectric actuator patches. Using the Heaviside function 
that corresponding to the inplane location of patches, the strain induced by the actuators 
was represented. In addition, the equivalent external forces induced by the piezoceramic 
patches under some voltage field can be determined upon the assumption of free constraint 
for expansion or contraction of the patches. 

In the work of Heyliger and Saravanos [59], discrete- layer theories were developed 
to analyze the laminated composite beams and plates which contain active piezoelectric 
layers or patches. The coupled relationship between the elastic and electric variables were 
explicitly represented in the governing equations. Both the static and dynamic behaviors 
were also considered. Later, Heyliger and Brooks [60] derived the exact solution for 
piezoelectric laminates in the two-dimensional configuration of cylindrical bending. These 
results provided useful information to evaluate the developed or future piezoelectric plate 
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theories. 

Development of piezoelectric shell theories were initiated by a number of Russian 
authors, whose efforst are summarized in the outstanding monograph of Kudrayavstev 
and Parton. Yet only recently have these been extended for active control problems. 

In 1989, Tzou et al. presented a numbers of works regarding piezoelectric shell theory 
[42, 49, 51, 55]. In [42], a laminated thin shell with piezoelectric layers was evaluated. 
The governing equations of the dynamic state were derived based on Love’s hypothesis - 
and Hamilton’s principle. Later, Tzou [49] introduced a piezoelectric element to analyze 
the distributed sensing and active vibration control of flexible plates and shells. The finite 
element formulation was also given. However, only a zero-curvature shell problem was 
demonstrated. In [51], Tzou and Tseng developed a “thin” piezoelectric solid element with 
the internal degrees of freedom. This element was used in the finite element formulations 
to analyze the piezoelectric shell. In 1991, Tzou and Zhong [55] developed electrornechanic 
equations of motion of generic piezoelectric shell using Hamilton’s principle and linear 
piezoelectricity. Numerical results using finite difference technique were represented, as 
well as a comparison with the experimental results. 

In addition to the works of Tzou et al., Lammering [54] developed a shell deformation 
finite element to analyze a shell structure with surface coated piezoelectric layers. In this 
work, a shear deformation elastic shell theory of the Reissner-Mindlin type was used to 
develop the finite element formulation. 

The Linear Theory of Piezoelectricity 

To begin the analysis for laminated structures with embedded piezoelectric layers, the 
material properties and governing equations for piezoelectric media must to be defined. 
The behavior of piezoelectric material used in the composite is assumed to be linear. A 
brief introduction to the material behavior and related equations is given below. 

Behavior of Piezoelectric Material 

In 1880, the Curie brothers discovered the piezoelectric effect on certain crystalline 
materials. When mechanical stresses are applied on a dielectric crystal, surface electric 
charges are instantaneously created. The converse effect is that the presence of an electric 
field results in changing the shape of the crystal. More than a hundred ferroelectric ma- 
terials with the piezoelectric effect have been found. Conventional piezoelectric materials 
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are crystalline or polycrystalline style. Therefore, they are mostly produced in the form 
of ceramics or crystals. The brittle properity leads to disadvantages in usage. 

In 1969, a strong piezoelectric effect in PVDF, polyvinylidene fluoride polymer, was 
discovered by Kawai [84, 16] . PVDF-type materials are flexible. The polymeric properi- 
ties are very different to the properities of the conventional piezoelectric materials. They 
can be produced and formed into thin-film sheets that are easy to cut or shape for com- 
plex configurations. Because of the convenience, structures with highly distributed sen- 
sors/actuators (intelligent structure) can be constructed. 

In laminated structures, the embedded or coated piezoelectric layers function as dis- 
tributed sensors and actuators. When applied forces result in strains within the structure, 
the surface charges on each layer can be collected through a surface electrode (i.e. it can 
be coated on the layer) to an outside detector. The deformation can therefore be mea- 
sured. Conversely, if a voltage field is linked to a certain piezoelect ric layer, the shape 
change within that layer can be used to actuate the structure. The surface electrode, 
which covers the piezoelectric layer, can be placed at any desired location. Only the 
portion of the piezoelectric layer covered by the electrode can initiate the effect. 

In Figure 1, a piezoelectric film with covered electrodes, and the poling and rolling 
(stretching) directions are shown. Figure 2 presents three possible deformations of piezo- 
electric laminated plate when applying an electric field. Both figures are taken from the 
work of Lee [84]. 

Equations to Describe the Piezoelectric Medium 

To describe the behavior of a piezoelectric medium, there are five mechanical and 
electromagnetic relations (or equations) involved : 

(1) the stress equations of motion/equilibrium, 

(2) the strain - displacement relations, 

(3) the charge equation of electrostatics, 

(4) the electric field - electric potential relations, and 

(5) the constitutive equations. 

These are given sequentially as follows: 

Tij,i = pU j (fj = 0) (1) 

Here r,j are stress components, p is material density, Uj are components of acceleration, 
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and the body forces fj are assumed to be zero. 

Sij = 2 ( u «.i 4" U J>>) 

Here Sij are strain components and U{ are components of displacement. 

Di,i = 0 


( 2 ) 

(3) 


Here D{ are the components of the electric displacement, or electric flux density, along 
the i direction. 

Ek — —<P,k (4) 

Here Ek are electric field along k direction and (f> is electrostatic potential. 

The piezoelectric constitutive equations present the relations among : 

(1) the stress field, the strain field, and the electric field; 

(2) the electric displacement, the strain field, and the electric field. 

They can be written as follows. 


Tij — CijklSkl &kijEk 


(5) 


Ei — e iklSkl T kEk 


( 6 ) 


The elements Cijki are elastic constants, ekij are piezoelectric coefficients, and £,•* are 
dielectric constants. In this research, the constitutive relations of an orthotropic layer off- 


axis are followed. The elastic stiffness matrix and the matrix of piezoelectric coefficients 


are given as 


and 


Cn 

C 12 

C 13 

0 

0 

C16 

0 12 

C 22 

C 23 

0 

0 

C26 

C 13 

C 23 

C 33 

0 

0 

C36 

0 

0 

0 

C44 

C45 

0 

0 

0 

0 

C 54 

Css 

0 

Cl6 

C26 

C36 

0 

0 

C66 


' 0 

0 

0 

ei4 

^15 

0 ‘ 

0 

0 

0 

e24 

^25 

0 

. e 31 

e32 

633 

0 

0 

^36 


(7) 


( 8 ) 


The indicial notations C^ki change into C a j, also e^j change into e* u (a, 6 = 1 , 2, • - - , 6, 
and k = 1,2,3). 
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The total set of 22 equations in 22 variables can fully define a given piezoelectric 
medium. Several additional equations related to linear piezoelectricity are briefly intro- 
duced here. 

The equation to describe the total surface charges q(t) at time t, which are collected 
by a surface electrode and sensed by the outside detector, is given by 


q(t ) = J D • ndS = J D^dS 


( 9 ) 


Here n is the normal vector of the small surface area dS, and D is the electric displacement 
vector in the piezoelectric medium. This equation represents that the D passing through 
the surface with electrode are collected. The internal energy function can be written as 


V = J y {\cv,SijS u + \e<iE,Ei)dV 


( 10 ) 


The first item of right hand side is the strain energy per unit volume. The second is 
electrostatic energy density. Finally, the electric enthalpy H per unit volume is defined 


as 


H{ S,E) = U-ED = U -EiDi 

— ^ CijklSij Ski &kij E k S v h/ t bj 

The electric enthalpy will be used in Hamilton’s principle in the next section. 


(11) 


Hamilton’s Principle for Linear Piezoelectricity 

In this section, the generalized Hamilton’s principle for a piezoelectric medium is pre- 
sented. The variational formulations developed in subsequent chapters depend critically 
on this expression. For the plate geometry, this expression will be given in terms of rect- 
angular Cartesian coordinates. For the shell, the variational formulation will be developed 
in both the cylindrical and the curvilinear coordinates. 

In this research, the generalized Hamilton’s principle is used as 

rt! rt 1 


Where 


[' 6(fC-U)dt+ [ ' 8Wdt = 0 
Jto Jto 

K.ffa-w 

6W = J v (f Su + <5(E • D ))dV + J { t <5u - a6<j>)<lS 


( 12 ) 

(13) 

(14) 
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Here K. is the kinetic energy of the system, U is the internal energy of the piezoelectric 
medium, and <5W, virtual work done to the system. The four items on the right hand 
side of equation 14 represent the virtual work done by (1) body force in a virtual dis- 
placement, (2) variation of electric field in a electric displacement, as- well as variation 
of electric displacement in a electric field, (3) the prescribed surface tractions in a vir- 
tual displacement on medium surface, and (4) the prescribed surface charge <7 in a small 
variation of electrical potential <f>. 

Since the applied electric field E and electric displacement D are interdependent, the 
system is not conservative. The body force f is assumed to be zero here. The expression 
of equation 14 is the governing equation used in this research, and is written as 


6 [ dt [ \]-puiUi — H(Ski,Ek) dV + [ dt [ (ti8ui — cr8<?)dS = 0 (15) 

Jto jv L2 Jto Js 


Here H(Ski,Ek ) is the electric enthalpy, which replaces the terms of internal energy and 
E • D in equation 14. The value will not be changed due to the coordinate transformation. 
The variational form of H is expressed as 


8H = 8(U - E • D) = SU - 8{ E • D) 

= CijkiSij8Ski — ekijEkSSij — ekij8EkSij — e l] E x 8E ] (16) 

These equations will be used extensively in later chapters. Depending on the problem 
geometry and coordinate system used, the form of these equations will be different and 
are detailed for each case in the following chapters. 
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3. SEMI-ANALYTIC SOLUTIONS 


This section details a discussion of what are termed semi-analytic solutions for simply- 
supported laminated piezoelectric plates. This work was an extension of work completed 
by the principal investigator and the technical monitor of this research project (DAS). 
The phrase “semi-analytic” denotes the fact that use of Navier-type displacement fields 
will identically satisfy the in-plane requirements of a simply supported plate for the exact 
solution of the equations of linear piezoelectricity. In this case, as for the exact solution 
included in a later chapter, the problem effectively becomes one-dimensional. 

The new and major thrust of this section is to detail 1) the relative accuracy of the. 
two major out-of-plane displacement approximations used in previous work in comparison 
with exact solutions, and 2) description of an alternative displacement field that was 
found necessary to accurately model piezoelectric plates. This field uses piecewise linear 
variation of the out-of-plane displacement components that are different than those for the 
in-plane and potential components. This is discussed to some extent in the chapters on 
the plate and shell elements, but is included here because it was for this class of problem 
that a preliminary numerical algorithm has been completed and results computed. 
Governing Equations 

Geometry 

The geometrical configuration of the laminate is such that the thickness dimension of 
the laminate coincides with the z-direction, with the lengths of the plate in the x and 
y directions denoted as L x and L y , respectively. The general problem considered in this 
study is to determine the behavior of the elastic and electric field components throughout 
the laminate under an applied mechanical or electrical loading. Each layer of the laminate 
can be composed of a purely elastic, piezoelectric, or conducting material. The forcing 
function is introduced through either an applied surface displacement, traction, potential, 
or electric charge. 

Variational Formulation 

A single piezoelectric layer has the constitutive equations given by [67] 

crij = CijSj - ejiEj ( 17 ) 

D{ — eij Sj tij Ej 
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Here <j,j are the components of the stress tensor, C tJ are the elastic stiffness components, 
Sj are the components of infinitesimal strain, eji are the piezoelectric coefficients, E, are 
the components of the electric field, D; are the components of the electric displacement, 
and are the dielectric constants. The poling direction in this study is coincident with 
the x 3 or z axis. 

The strain-displacement relations are given by 


c _ 1 ( duj , 8uj' 

~ n l + 


(18) 


2 \dxj dxi J 

Here S tJ are the components of the infinitesimal stain tensor and u, represent the displace- 
ment components. To be consistent with equation 17, the conventional notation for the 
strain indices has been used, i.e. Sn = Si, S23 = S4, etc. The electric field components 
E t - are related to the electrostatic potential <j> using the relation 


Ei = - 


d<t>_ 

dx{ 


(19) 


For the materials used in this study, it is assumed that the non-zero components of the 
rotated piezoelectric tensor e,j are e 3 i, e 32 , e^, e 2 4, eis, e 2 s, eu, and e^. The elastic 
stiffnesses C,j are those of an orthotropic material rotated about the z axis, and the 
dielectric constants are given by Cn, ei 2 , e 22 , and e 33 . 

The starting point for the variational formulation is Hamilton’s principle for a piezo- 
electric medium [67], expressed as 


: />/J 


pujUj - H(S k i, E k ) 


dV + 



8uk — cr6g)dS = 0 


( 20 ) 


Here t is time, V and S are the volume and surface occupied by and bounding the solid, 
t and cr are the specified surface tractions and surface charge, respectively, 8 is the vari- 
ational operator, the . superscript represents differentiation with respect to time, and H 
represents the electric enthalpy. The electric enthalpy is given by 


H = ^CijkiSijSki - eijkEiSjk - \tijEiE } 


( 21 ) 
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The weak form of the governing equations, as well as the governing differential equa- 
tions themselves, can be found by applying the variational operator in equation 20 over 
a typical element. For the material constants of a typical lamina used in this study, the 
first variation of the electric enthalpy is 


8H = C 


£\$8E\ 




' d8v 

dz 


C 
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du d8u 

n di!h 

dv d8u 

' n fylh 

dw d8u 
dz dx 


+ C 
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du d8v 
dx dy 
dv d8v 


+ C 


13 


'22 


+ c ; 


23 


dy dy 
dw d8v 
dz dv 


+ C 


23 


+ C- 


33 


du d8w 
dx dz 
dv dSw 
dy dz 
dw d8w 

~dl~d7 


+ C 


26 


du . 

( d8u 

'~di 1 


dv , 

( d8u 

% 

[dy 

dw j 

( d8u 


+ 


d8v 


+ 


C. 


44 


-'45 


( dv dw 
[d~z + ^ 
' du dw 


u6v dSw 
dz + dy 
f d8v dSw 
dz dy 


dv ' dw 
dz + dy 


' 36 


' du dv 

, dy dx 




+ 


+ 


+ 


du dv ' 
dy + dx. 



- ei S 8Ei — 



d8w\ 


fr , f dv dw\ (du diu\ 

- e M SE 2 + e, s 6E, + g x ) - 

du dv ()w (du dv\ 

~ c * SE *Ty - e * iLi T; - x ' u [aZj + dx ) ■ 
t\\E\8E\ — C22E28E2 — t^E^SE-i — t\2Ei8 E2 — t\2E28 E\ — 

d8w\ 


( dSu d8w\ ( d8v d8w\ 

e ' sE ' (' IT + ir) " e ' ,E * ( a7 + WJ 


_ , d8u 

— e2s^2 I h 


dz dx 


d8u d8v d8w 

631 1 - e 32 e33 ^ 3_ ^Z e 36-C/3 


d8u d8v 
dy dx 


( 22 ) 


Using the assumption of periodic motion, the substitution of this expression into equation 
20 yields the final weak form, which provides the basis for the finite-element approxima- 
tions over an element. It is possible to integrate this expression by parts to give the three 
equilibrium equations and the conservation of charge equation. For brevity, results of this 
step are not included here. 


Displacement Functions 


There are several cases in both static behavior and free- vibration analysis (or which 
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either exact solutions are available or global in-plane functions can represent the elastic 
and electric fields within the laminate. One such case is for simply-supported plates. In 
this case, the Navier-type boundary conditions along the edges of the plate are given as 


w(x,ti,z) = w(x,L y ,z ) = w(0,y,z) = w(L x ,y,z) = 0 (23) 

<£(x,0,z) = <f>(x,L y ,z) = <£(0 ,y,z) = <p(L x ,y,z) = 0 

u(x, 0,z) = u(x,L y ,z) = 0 
v(0,s!,*) = v(L,,y,z) - 0 

(24) 

These conditions are satisfied exactly for the functions 

u(x,y,z ) = U(z) cos px sin qy = U exp(sz) cos px sin qy (25) 

v(x,y,z) = V(z)sinpxcosgp = V exp (sz) sin px cos ^v/ 
w(x, y,z) = W(z) sin px sin qy = W exp(sz) sin px sin qy 
(f>(x,y,z) = $(z) sinpxsin^y = $ exp(sz) sin px sin qy 

Finite Element Approximation 

The through-thickness approximation for each displacement component and the elec- 
trostatic potential can be given as 

U (M) = 

i=i 

»(*.*) = EW)*;w 

i=i 
i= i 
i=i 

(26) 

Hence the variables are represented by linear combinations of a known through-thickness 
distribution described using the one-dimensional Lagrangian interpolation polynomials 
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^j(z) [14,16]. The in-plane distributions are already included in the total approximation. 
Here (n-1) is the number of subdivisions through the laminate thickness (typically taken 
equal to or greater than the number of layers in the laminate), and U Jt - is the value of 
component u at height j corresponding to the i-th in-plane approximation function. 

The difference between the two theories considered here is in the form of 'I'". For 
the case of constant transverse displacement through the thickness, this function is equal 
to 1. For a variable-w approximation, this function can be any Lagrangian interpolation 
polynomial. It is also possible to use functions in the z-direction which are non-zero only 
over specified regions. This would be most useful for the functions for conducting 
materials, in which case the potential is a constant. Also, independent approximations 
for all variables can be used by denoting different ranges and limits of integration for each 
component. 

Substituting these approximations into the weak form and collecting the coefficients 
allows the governing equations to be expressed in matrix form as 


[M u ] [ 0 ] 

[0] [M 22 ] 

[ 0 ] [ 0 ] 

[ 0 ] [ 0 ] 


[0] [0] 1 
[ 0 ] [ 0 ] 
[M 33 ] [0] 

[0] [0] J 





+ 


[!<"} 

[A' 21 ] 

[A" 31 ] 

un 


[K 12 ] [A’ 13 ] 
[A 22 ] [A 23 ] 
[A 32 ] [A 33 ] 
[A 42 ] [A 43 ] 


[A**] 

[A 24 ] 

un 


1 

{0} 1 


) 

( y } 

► < 


1 {M / ) 


J 1 

i w J 



{A 1 } 

{A 2 } 

{A 3 } 

{<?} 


The elements of these matrices contain additional submatrices whose elements are de- 
termined by evaluating the pre-integrated elastic stiffnesses, piezoelectric coefficient, or 
dielectric constant through the thickness multiplied by the various shape functions or 
their derivatives as determined by the variational statement. These are identical to those 
of the plate element, and as such are not stated here. They are documented in the next 
chapter. 

The nature of the submatrices depends on the approximation used for w. For variable 
w, the structure of [K-' 3 ] are similar to those of the other matrices. Depending on the 
approximation functions used, this type of model is either similar or identical to a three- 
dimensional finite element model of the equations of piezoelectricity. This aspect and other 
computational issues regarding these types of models for elastic laminates are discussed in 
[61]. For constant w, the submatrices within [K 13 ], [K 23 ], and [I\ 43 ] are column vectors and 
those in [K 33 ] become scalars. In general, the submatrices within [K‘ J ] are each of order 
(n+1), while the [K u ] themselves depend on the order of in-plane approximation. These 
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can be defined by the subscripts a and /?. The final representation of these elements can 
be expressed in fairly compact form and are given in the chapter on the plate element. In 
matrix form, they are written as 


[M] [0] 

[ 0 ] [ 0 ] 




’ ftfAA 



'k** 





{ 0 } \ 
{ 0 } / 


(28) 


Assuming periodic motion and eliminating the potentials results in the eigenvalue 
problem 


([j?]-u 2 [M]){A} = {0} 

where 


[K] = (29) 

The numerical results for the preceding theories are given in a later chapter. However, 
it is necessary to note at this point that the case of a constant transverse displacement 
gave results for both the static and dynamic cases that was wholly inadequate for most 
geometries. Part of the reason for this was the fact that the transverse strain must be zero 
in this case, which implies that the actuation strain is also zero. For any dimension plate, 
this is a hindrance far more severe than for the pure elastic case, as this strain figures 
significantly in the stress computation and in the evaluation of the electric enthalphy 
terms. 

This type of behavior can be asauaged by introducing approximations for w that are 
at least linear through the thickness of the entire laminate. This is accomplished using 
an approximation function for w through the thickness of the laminate that is completely 
independent of the two in-plane components and the electrostatic potential. 

This type of approximation can be thought of in the following fashion. Each indi- 
vidual layer can be described by one or more “real” layers for purposes of finite element 
discretization through the lamina thickness, composed in this study of linear Lagrangian 
interpolation polynomials. Spanning over these layers is one or more "pseudo-layers” (see 
Figure 3), which is the range of layers encompassed by the linear approximation for any 
of the variables, which in this case will be the transverse displacement w. For example, 
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consider a plate with three physical layers, which is modeled using three layers (and 4 
nodes) for u, v, and ^ but a single pseudo-layer for w. Each of the linear basis functions 
for u, v, and <j) has support only over two layers, with the basis function for w spanning 
over all layers. 

The shape functions for the i-th layer for each of the true layers are given in terms of 
the global z-coordinate as 


^ = L h~ (30) 

Here h,- is the thickness of each of the true layers, and the t superscript indicates “true” 
and corresponds to the approximations for u, v, and <f> (and any, all. or none of these may 
have this type of approximation). 

The approximation functions for the transverse displacement component are given by 


= fit ~ oc— 

Lk 

ipl = 7i + Q 'J^ ( 31 ) 

Here is the thickness of the k-th pseudo-layer, and the values for a, 7 , and change 
for each true layer, with their values being given by 



(32) 


7* = Q «' ( 33 ) 

i=i 


A = l-i>i (34) 

J = 1 

By analytically integrating the products of the appropriate shape functions and their 
derivatives, the appropriate sub-matrices similar to those computed earlier for the original 
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two theories can be obtained. This approach has the significant advantage of allowing a 
transverse normal strain through the laminate thickness but introducing only (potentially) 
one new degree of freedom. By allowing a transverse normal strain, the difficulties of the 
constant-w theory disappear. Examples of this type of behavior are given in a later 
chapter, but indicate one of the major conclusions of this study: it appears that the 
constant-w theory is wholly inadequate for modeling the behavior of piezoelectric solids 
because of the inability to model the normal actuation strain. 
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4. THE PLATE ELEMENT 


In this chapter, the theoretical formulation for piezoelectrically laminated composite 
plates are presented. The variational form is obtained using Hamilton’s principle for a 
linear piezoelectric solid [67]. The equations resulting from this variational statement will 
be solved using the discrete-layer theory similar to that described in the preceding chapter, 
in which the three-dimensional elasticity theory is reduced to a two-dimensional laminate 
theory by assuming an approximation of the displacements through the thickness. 

In this study, a piece-wise linear variation is assumed for the in-plane displacement and 
for the electrostatic potential components through the thickness. Regarding the trans- 
verse displacement components, three theories will be presented: 1 ) a constant transverse 
displacement through the thickness, 2) a piece-wise linear variation through the thick- 
ness, and 3) an independent piece-wise linear variation through the thickness different 
than that of the in-plane and potential components. 

Governing Equations 

Electric Field- Potential and Strain- Displacement Relations 

The electric field components are related to the electrostatic potential by [67] : 


Ek = — v?,jt 

In Cartesian coordinates, the above expression yields: 


(35) 


F - F - 

El dx E 2 dy 


E 3 = - 


dp 

Ik 


(36) 


In addition, the strain-displacement relations are given by: 

1 


(37) 


Stj ~ cy ( U lJ "t" U j,‘) 

which turns into the following six components when is expressed in cartesian coordinates: 

r 

:) 


du 

dx 

1 du dv, 
xy ~ 2 [ dy + dx‘ 

dv 

1 dv dw, 
yz ~ 2 [ dz + dy 

dy • 

dw 

„ i .du dw 

S - = 2<fc + fc 

~k 
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Constitutive Relations for a Piezoelectric Material 

For a piezoelectric material, the mechanical behavior is coupled together with the 
electrical behavior. This electromechanical coupling is described by [67]: 

= CijkiSu — CkijEk (38) 

D{ = eu-iSki + €« kEk (39) 

Here <7,y represent the components of the stress tensor, e^y are the components of the 
dielectric permittivity or the piezoelectric constants of the solid, e^. are the components 
of the dielectric tensor, D{ are the components of the electric displacement tensor, and 
Cijki are the components of the elastic stiffness tensor. 

Since we are interested in a linear theory, it is important to note that: 

Cijki — Cijlk — Cjikl — C/ klij 
&ijk = & ikj 
e ij = e ji 

Taking into account this symmetry of the material, Nye [68] compressed this notation 
by replacing ij or kl by p or q, where i,j, k, and l take the values 1,2. and 3, and p and q 
take the values 1,2, 3, 4, 5, and 6. Hence the equations can be rewritten as: 


CpqSq ek p Ek 

p = 1,2, — ,6 

q= 1,2,-.. ,6 

(40) 

CiqSq “h 6ikEk 

i = 1,2,3 

k = 1,2,3 

(41) 


Now that the elastic, piezoelectric, and dielectric constants are specified by two indices, 
they can be written in matrix form as follows: 


Cu 

C U 

C 13 

C 14 

C\5 

Cie 

C\ 2 

Cyi 

C 23 

C 24 

C 25 

C 26 

C 13 

Crs 

C 33 

C 34 

C 35 

C 36 

C 14 

C 24 

C 34 

C 44 

C 45 

C 46 

C 15 

C 25 

C 35 

C 45 

C55 

C$ 6 

. Ci 6 

C 26 

C 36 

C 46 

C 56 

C 66 



^11 

ei2 

613 

614 

^15 

6 16 

Cfcp = 

^21 

622 

623 

624 

^25 

626 


^31 

632 

633 

634 

£35 

636 


(42) 


( 43 ) 
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C ij “ 


(44) 


£ 11 e 12 e 13 

£ 12 £ 22 £ 23 

£ 13 £ 23 633 

At this point, we have 21 independent elastic constants, IS independent piezoelectric 
constants, and 6 independent dielectric constants. Because of the symmetric properties 
of some piezoelectric materials, these matrices can be further simplified. For instance, for 
a monoclinic material, which can represent an orthotropic layer oriented off-axis, these 
matrices can be reduced to: 



On 

C12 

C13 

0 

0 

C \& 


Ci 2 

C22 

C 23 

0 

0 

C 26 

C 13 

C 23 

C 33 

0 

0 

C36 

0 

0 

0 

C44 

C45 

0 


0 

0 

0 

C 45 

C 55 

0 


Cie 

C26 

C36 

0 

0 

C 66 


0 

0 

0 

ei4 

615 

0 


0 

0 

0 

624 

625 

0 


. e 31 

£ 32 

633 

0 

0 

636 



(45) 


(46) 


£ 




e u 0 0 

0 £22 0 

0 0 £33 


(47) 


Thus the material with this type of symmetry is described by 13 independent elastic 
constants, 8 independent piezoelectric constants, and 2 dielectric constants. 


Variational Form of the Governing Equations 

In classical mechanics, the Hamilton’s principle for a system with nonconservative 
forces states that: 

6 [ l Ldt + /* 8Wdt = 0 (48) 

Jto Jto 

where 

L = r(x fc ) - V(xk,t) = L(ik,Xk,t) 

In the above expressions, t is time, L is the Lagrangian energy function, T is the kinetic 
energy, V is the potential energy, 8W is the work done by the noncoiiservative forces in a 
virtual displacement, and 8 is the variational operator. 
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As noted in earlier chapters, Hamilton’s principle for a zero body and surface forces 
is given by 

s 1 dt J ^ /5UjU ' i “ - eijEiSj - ^ CijEiEjldv = 0 (49) 

The weak form of the governing equations can be found by substituting the elastic strain- 
displacement and the electric field-electric potential relations and by operating the vari- 


ational operator 8 on Hamilton’s principle. Assuming static conditions, the above proce- 
dure yields the following for an off-axis, orthotropic piezoelectric layer: 


du dSu du d8v du d8w du d8u du d8v 

on = Cut ; — 5 |- + W 3 — — — b ^i 6 ("X — x b zz — t — ) 


dx dx 


' dx dy 


'dx dz 


'dx dy dx dx 


n dv d6u ( n dv dSv n dvdSw n ( dvObu dv dbv v 

b C 22 T . — ~ b C 237 ; — x b Cm! ~ — r: b ~ — x — ) 

dy dx dy dy dy dz dy dy dy dx 


dz dx 


dz dy 


dz dz 


dz dy 


dv dSv dv dSw dw d8v 

44 dz dz + dz dy + dy dz 

dv d8w dv dSu dw dbw 

45 dz dx + dz dz + dy dx + 

dw d8v dw d6w du dSv 
+ ^ + — — + 


C 45 ( 


dx dz dx dy dz dz 


dw dSw dw dSu dudSw 

^55 \ o r\ + n n + A A *b 

dx dx dx dz dz dx 
_ .dudSu dvdSu dudSv dvd8v s _ .dudSw 

+ stst) + Cxls-s- + t-^-) + c x (—— + 


dy dx dx dx 


dy dy dx dy 


dy dz 


du d8u du dSv dv d6u 
66 dy dy dy dx dx dy 

dv dSip dw dSip . . dw d8p> 

e »( AI-3 T + at ^r) + e. 5 (— ^ + 


dz dx 

dw dSw 
dy dy 

diu dSu 
dy dz 

du dSw 
dz dy 

du dSu 
I~z~dl 

dv dSw 
dx dz 

dv dSv 


'dz dx dy dx 


dx dx 


. dv dSip dw dSip dw dSip 

+ a—a^) + e 2s(“ ^ + 


1 dz dy 


du d8ip dv dSip 

e 3i ~ ~ — b 632 - 5 — ^ — b 633 


dy dy 
dw dSip 


dx dz 


dy dz 


dz dz 


dx dy 

du d8<p 

+ e 36 ( — 7T~ + 


dx dx 

du dSip 
dz dx 

du d8<p 
dz dy 

dv dSip 


dy dz dx dz 


. dip d8v dip d8w . , dip dSw d'p dSu 

+ Sisr) + e ls (7r-— + 


dx dz dx dy 


dx dx dx dz 


+ 


+ 


_ dw dSu dw dSv dw dSw . dw d8u dw d8v . 

C xz — — + C 23 — — + ^33 - — — + Cm( — — + -t-tt) + 


) + 


) + 


) + 


) + 


) + 


) + 


) + 


) + 


) + 


) + 
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. dp dSv dp dSw dp dSw dp d8u 

«»(£— + ■&—) + Mrx- + «r-sr) + 


' dy dz 


dp dSu dp dSv 
631*7; — 7 1 - 632-7; — 7 — h £33 


dy dy 
dp dSw 


dz dx 


dz dy 


dz dz 


1 dy dx dy d 

. dp dSu dp d8v 

+ 636 V 7 b 


dz dy 


611 


dp d8p dp dSp 


dx dx 


— C22 


dy dy 


~ «33 


dz dx 

dp dSp 

Iklk 


) - 


( 50 ) 


Discrete-Layer Approximation 

Three theories are presented in this work. The first asssumes a piece-wise linear vari- 
ation of the transverse displacement through the thickness. The second yields a constant 
transverse displacement through the thickness. The last uses independent piece-wise lin- 
ear approximations for w and the remaining in-plane and potential components. All 
theories are based on the general laminate theory of Reddy [32] for elastic laminates, with 
the added feature that the electric potential is included as an additional variable. 

A piece-wise transverse displacement model 

This theory is based on approximations of the displacement and potential variables in 
the following form: 

n m n 

u(x,y,z,t) = «)*>(*) = EE £««)«•(*. »)**(•=> 

j= 1 *= 1 j = 1 

n m n 

»(*.» a,0 = £ = ££ 

j = 1 *= 1 i=i 

n m n 

w(x,y,z,t) = '£W i (x,y,t)^(z) = »)*■(*> (51) 

j=l «= 1 j=l 

In similar way, the approximation for the potential can be written its: 

p(x,y,z,t) = <)$?(*) = ( 52 ) 

j=t «=i i=i 

where u,u, and in represent the displacement components in the x,y , and z directions 
respectively of a material in the underformed laminate, and p represents the electric 
potential. In the above equations, two approximations have been made. In the first one, 
the transverse variation of the displacement field is defined in terms of the one-dimensional 
Lagrangian interpolation polynomials 'Pj(z) which are associated with the j th interface 
of the layers through the laminate thickness and are defined only on two adjacent layers. 
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The laminate is divided into n nodes distributed through the thickness at the j th level, 
defined by z = Zj in the undeformed laminate. At this level of approximation, the 
functions Uj(x,y,t),Vj(x,y,t),Wj(x,y,t) , and $j(x,y,t) represent the displacement and 
potential components of all points located at the j th plane. In the second approximation, 
the in-plane displacement and electric potential field have been defined in terms of the 
two-dimensional Lagrangian interpolation polynomials ^(x,y). The j i>l plane is divided 
into two-dimensional finite elements and m represents the number of nodes per every one 
of these elements. At this level of approximation, the functions Uji, Vj,-, Wji, and 
represent the values of the components u,u,and w , and <p at height j corresponding to 
the i th node of the two-dimensional finite element. 

Now substituting equation 51 and 52 into the variational form of Hamilton’s principle 
50, integrating with respect to the thickness coordinate z , and collecting terms, the final 
equation in matrix form can be expressed as: 


[*].{*}. = {F}< 


(53) 


or 


‘ [I< n ] [A 12 ] [A 13 ] [A 14 ] - 


{«} 1 


{/'} I 

[a 12 ] [ a 22 ] [a 23 ] [a 24 ] 


W 

i — j 

in 

[A 13 ] [A 23 ] [A 33 ] [A 34 ] 

\ 

M 

* 1 

in 

[A 14 ] [A 24 ] [A 34 ] [A 44 ] 

e 

{*} J 

e 

{/“} J 


(54) 


The elements of the stiffness matrix [K] e and force vector {F} e are given below: 


[K n U = [ 

J A 




dx dx 


/3\D U AvD u 

[D 55 ]^^ + [ A ™\ Qy - dy \ 


dy dx 
dxdy 


dx dy 


)+ 


I *”i- - L 


ia u ] 0*I? , 5 


(b«], + 1^«] 


+ 

dx dy 1 1 dy dy 

dVl !>% 




dx dx 


dy dx 

dm 


dxdy 


r T rl lT / u r) lT/ u f) lb ,v 

[I< 13 U = J a [ b 13 }-q^^ + [ 536 ]-^^ + [B 45 }'H“-^ 
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[K U U 


[K”U 


[I< 23 U 


[K 24 U 




[K*U 


[I<"U 


I B*)K-j£ 


dxdy 


-L 


[£“]«; ^ + [£ 1S 1 *:^ + [£ 31 ]^^+ 


dx 


dx 


r)<h u 


dxdy 


| [A 22 J^^ + [A 26 ](^-^ + 


dy dy 


dx dy dy dx 


7 ^)+ 


^vD v < 9 \I/” 
+ [>*“] _£ 


dxdy 


-L 


„ d<h v „„ /M" -_ r>V'V 

[B 23 ]^'!'? + [B 36 ] — -±m + [B 44 ]^* — - + 
dy p dx p dy 


[B 45 tf 


c W 


XV 


45 mu 

“ dx 


dxdy 


-L 


dM 


+ [^n-zr + [E 32 }^n+ 


dm 


32 


dm 


° dx 


dy 


dy 


fl\ T>u 


dxdy 


-L 


[A ^_n^n + [a « ]( ^5 + ^^i )+ 


dy dy 


dx dy 5?/ dx 


ftyhxv f)\h w 


dxdy 






d_m 

dx dy 


+ [B 33 ]^^ 


dxdy 


-L 


[G ] 


d^i d^p 
dx dx 


+ [G 22 ] 


dy dy 


+ [G 33 ]'P;/P* 


dxdy (55) 
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toti 


{<z} = 

/ D^idxdy 

Ja 

(56) 

{/“} = 

f fi^idxdy 
Ja 

(57) 


The matrices involved in the above equations result from the evaluation of the integra- 
tion of the elastic and dielectric constants through the thickness of the laminate. These 
matrices can be written as follows: 



(58) 

(59) 

(60) 

(61) 

(62) 

(63) 

(64) 


'-pkm 





Gkm 


d^jz) d^jz) 
dz dz 


dz 



£ km 


dVj(z) d^jjz) 
dz dz 


dz 


J=1 JZ ‘ 


(65) 

( 66 ) 

(67) 


where N is the number of subdivisions through the thickness and, in general, it will be 
taken equal or greater than the number of layers in the laminate. 

The theory which allows independent piece-wise linear approximations for w and the 
remaining components can be represented by the above discussion, with the provision that 
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now the through-thickness linear Lagrangian interpolation polynomials will different. 

A constant transverse displacement model 

The following approximations for the the displacement and potential variables are 
assumed: 

u(x,y,z,t) = Y / U j {x,y,t)^(z) = £ £ 

i=l t=l j = 1 

v(x,y,z,t) = 

i=i «=i j=x 

*=1 

n _ m n 

<p{x,y,z,t) = (a-, j/pl'/U) (6S) 

j=i «=i j = 1 

The only difference between these approximations and those assumed in the other models 
is found in the approximation used for the transverse displacement in the r direction. In 
this model, the transverse displacement is considered constant through the thickness and 
therefore the one-dimensional Lagrangian interpolation polynomials ^(s) used before is 
equal to one. The represents the transverse displacement corresponding to the i — th 
node of the two-dimensional finite element used in the approximations through the x — y 
plane. 

Proceeding in the same way as that used in the other model, the following element 
stiffness submatrices of \K\ e and the force vector {jp} e are obtained and may be written 
as: 


l* U l 


a(3 


- L 


[A 11 ] 


dx dx 


+1 ^ l( ^s + M 

dy ox dx dy 


)+ 


<9# o 


dxdy 


[I< U U = / 

J A 


1,4 1 dx dy 


fix i» u d 'I/o c/vi i" 

+ [a 26 ]^^ + [a 16 i v 4t^ 

a - ox dx 


dy dy 


+ 


fixhu d^ V o 


dxdy 
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[*“]«* 

[I< 14 U 

[K 22 U 


[K 2Z U 

[I< 24 U 

[I^U 

[K M U 

[K 44 U 




0 , { A D S5^ lTl u9^W0 


^\ {AB « )n ^ + {AB s> )K & 


dxdy 


i.[ 


[£ 25 ]*“^2 + [£ 1S ]^^2 + [£ 31 ]^ f ^+ 


dx 


dx 


fl\ Ir“ 


dxdy 


L 


ahOTSWj 


[A 22 ] 


+ [A 26 ]( 


dy dy 


,d$”d^ d^d^ 


-| 2 . — 

dx dy dy dx ' 


F)^ v d^\ 

[D ^ W + [A eef_^_J 


dxdy 


L 


{AB 44 }V v q ^- + { AB 45 }^ v a dK]!W 0 


d^ 

dy 


dx 


dxdy 






dxdy 


UA «™i d Ji+ AA « { 


dy dy 


dV%dV% d^d^ 


dx dy dy dx 


AA 


„d®jd*j 


dx dx 


dxdy 


\(EE ) 


d^d% 15 d*«d*£ w d^d^ 


dy dy 


+ [EE 15 ) 


dx dx 


+ (££")• 


dy dy 


■ + 


(EE 25 ) 


dW? d^l 


dx dy 


dxdy 


if 


np? 

A I OX OX 


+ (G 


22 l d^d^ 

dy dy 


+ [G 33 ]*^ 


dxdy 


( 69 ) 


( 70 ) 
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As before, the matrices and vectors involved in the above equations result from the 
integration of the elastic and dielectric constants through the thickness of the laminate. 
These matrices and vectors can be expressed as : 


At 

N 

= E f +1 c fcm ^,(z)^(z)dz 

;=i Jz > 

(71) 

ot 

jrr*'c 

J Z[ dz dz 

(72) 

AB* m 

= 

J zi dz 

(73) 

E?r 

" p+l d$i(z) - . . 

= E / e ^m , *j[z)dz 

l~l J Z l aZ 

(74) 

Et 

" P+. - . .d*i(z) , 

— E / t d~ 

i = i Jz ‘ az 
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Transformation Matrices 




In this work, it is desirable to simulate the behavior of a shell element derived from the 
plate element using an assemblage of flat plates. In this sections it is assumed that the 
behavior of a continuosly curved surface can be adequately represented by the behavior 
of a surface built up of small plate elements. 

Consider a typical plate element as derived in this section of the report. In general, 
there are four degrees of freedom per node in the j th level, and depending on the approx- 
imation used for the transverse displacement, w, in the z direction we would have either 
4*(number of layers +1) when w is a linear function through the thickness of the laminate 
or 3*(number of layers +1) +1 when w is considered constant through the thickness of 
the laminate as the total number of degrees of freedom per node through the thickness. 
The independent approximation case can be considered in this section as well. 

The stiffness matrix derived before was based on a system of local coordinates, there- 
fore a transformation of coordinates from a local system, (x\y’,z ) to a global system 
(x,y,z) will be necessary in order to assemble the elements and to write the appropiate 
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equilibrium equations. Moreover, it will be convinient to read the element node coordi- 
nates in the global system and get from these the local coordinates. 

Therefore according to figure 14, we have that: 

{ V '} = pr]{V} (78) 


where V- are the components of a vector on the local system , Vi are the components of 
a vector on the global system, and Tij are the cosine angles between the local and global 
axes and the matrix is given by: 

" 1 0 0 ' 

0 cos/3 sin/3 (79) 

0 —sin/3 cos(3 

The cosine angles can be calculated by the following expressions: 

Vi ~ Vi 


T - — 

1 ij — 


cos/3 = 


sin(3 = 


-j 


(80) 


\/((*j - + (»j - ».■)’) \/((2j - *i ) 3 + (i/j - v.) 1 ) 

Since we have to transform both the forces and displacements by using the equation 
7S and taking into account that [&f oca /]{ A/ oca /} = {F; oca ;}, we get the final relationship: 

([rfWo^riHA,,,,.,} = (si) 

but 


= pf 

therefore we can write equation 81 as: 

{ ^global } {Fglobal} 


(82) 


(83) 


where A; are the components of the displacement vector, F, are the components of 
the force vector, k{j are the components of element stiffness matrix in local coordinates, 
and Kij are the components of the element stiffness matrix in the global coordinates. 

Therefore for one node in a linear plate element and for both theories, with w varying 
and w constant through the thickness, the corresponding displacement vectors will appear 
as 


U n } 
U n 
V\\ 
Vu 
Wn 
W n 
$11 

. $12 . 


( 84 ) 
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( 85 ) 


U n 

U 12 

Vu 

< V 12 1 

VTi 

$11 

. $12 „ 


and the corresponding transformation matrix [T] will be: 
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0 
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0 

0 

0 

0 

0 

0 

0 

0 

1 


Examples of applications of this element are included in the results chapter of this 
report. As of this writing, only the flat plate case had been investigated with the assembly 
of plates remaining to be studied. 
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5. THE SHELL ELEMENT 


In this chapter, the formulation for two laminated piezoelectric shell elements is given: 
the cylindrical shell and the general curvilinear shell. Both are based on solutions of the 
weak form of the equations of motion using the discrete-layer type theories described in 
the two preceding chapters. 

Variational Formulation in Cylindrical Coordinates 


A special case of a laminated shell is the cylindrical shell geometry, which is best 
formulated in cylindrical coordinates. Let u,v,w be the components of displacement along 
the r,0, and z directions in cylindrical coordinates, respectively. The strain-displacement 
relations and the stress equations of motion are given as 
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Because the normal vector of layer surface will be taken along the / -direction, 


also this 


direction is assumed to be the poling direction of each layer, the matrix of piezoelectric 
coefficient is given as 

fin ei2 ei3 0 0 ei6 


0 0 0 e24 ^25 0 

0 0 0 634 635 0 


(91) 


Also, the electric field can be written as 


E(r,d,z) 


dq> _ 1 8<j> _ 8$ _ 

or r dO dz 


( 92 ) 


Hamilton’s principle as writen in equation 14 can be expressed as 
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This state provides the basis for finite element approximations to the equations of 
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motion as formulated in cylindrical coordinates, and are a special case of that described 
below. This more general case is described thoroughly and has a corresponding numerical 
algorithm that was completed as part of this proposal. 

Variational Formulation in Curvilinear Coordinates 


To analyze the laminated piezoelectric shell with arbitrary geometry, curvilinear coor- 
dinates ( see Figures 4 and 11) will be used here to overcome the disadvantage from the 
integration directly using the global coordinate system. The local and global systems can 
be related using mapping technique. Let u, v, w be the components of displacement along 
the x-, y-, and z-coordinate directions of a Cartesian coordinate system, respectively. This 
coordinate represents the global coordinate system. The curvilinear coordinate system, 
described by £-, 7/-, and ^-coordinate directions, is the local system for the general shell. 
Let ti, v , ui be the components of displacement along the t/, ^'-coordinate directions. 
A laminated shell with complicated geometry can be separated into several small parts 
that function as shell elements, and described using more but simpler curvilinear coordi- 
nates. In this section, the variational formulation will be expressed in terms of the local 
coordinate system. 

The strain-displacement relations and the stress equations of motion using the global 
coordinate are given as 
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d<J Z y 

dy 
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(98) 


The electric field can be written as 


p( , d<f>_ d(t>^ d<f>_ 


(99) 


The strain components in global coordinate system can be changed into the local system 
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with the coordinate transformations given by 
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Equation 100 can be further expressed as : 
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The (u,v,u;) in x,y,z system and (u,v,w) in kViC system can be related using 
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The electric field (E x , E y , E z ) and (E^, E v , E zeta ) can be related using 
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The first differential of (u,v,w,<j))gi 0 b a i to (x,?/,z) and (u,v,w,(f>)ytob t i to (£,*/, C) can 

T 

related using the matrix, Jj -1 (f^)] • 


' u ^r 

u ,y 

U rZ 


’ U A 

u ,v 

n,c ’ 


§L 

dx 

*L 

dy 

1 

dz 

v ,x 

v ,y 

V ,z 


v ,c 

v,n 

v,C 


dr] 

dv 

dju 

W ,x 

W.y 

W ,Z 


W£ 


w ,c 


dx 

K 

dy 

K 

dz 

dc 

. fa r 

^»y 

fax \ 


. 4 >a 

t.T) 

fac \ 


_ dx 

dy 

dz . 


( 105 ) 


37 


The ^-direction is located along the normal direction (also the poling direction) of the 
layer surface, and the f- and 77- direction parallel to the layer surface. Hence equations 100 
and 102 can be used to present the elastic stiffness matrix and the matrix of piezoelectric 
coefficient in the curvilinear coordinate system. The variation of electric enthalpy, SH , 
can be expressed as the following equation. 
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Where 





[f/mj El] global = 

| ixx ^yy ^zz ^yz i-zx ^xy Ex Ey E 2 



(107) 

local = 1 

hi hv he hi hi hv E{ E n £ ( ] 



(108) 

<f>l] xyz = [ 

du du du dv dv dv dw dry du> — d<t> — d<t> 

dx dy dz dx dy dz dx dy dz dx dy 

-do 1 
dz \ 

(109) 

Kii&W = [ 

du du du dv dv dv dw dw dw —d<t> 

d£ dri d{ d£ drj d{ d£ drj dC, d £ 

— d<}) 

*T) 

-do 1 

J 

(110) 


■ J~'( g&) 0 0 0 

0 J-Hgf) 0 0 

0 0 J 

0 0 

1 0 0 0 0 0 0 0 

0 1 0 0 0 0 0 0 

0 0 1 0 0 0 0 0 

00020000 
00002000 
00000200 
0 0 0 0 0 0 1 0 

0 0 0 0 0 0 0 1 

00000000 


-V £2*1 
'ivi> 

0 

0 
0 
0 
0 
0 
0 
0 
0 

1 J 9x9 






ivi 1 -I 12 x 12 


(112) 



38 


[B} = 


100000000 

000001000 

000010000 

000001000 

010000000 

000100000 

000010000 

000100000 

001000000 

000000100 

000000010 

000000001 


( 113 ) 


12x9 


(H- 


[-^16x6 

0 


0 




/ J3x3 


in = [f •] [«] 


9x9 


[C] 


global 


= T 


[*] PL., [*) 


( 114 ) 

( 115 ) 

( 116 ) 


The matrix [ C]i oca i that include the elastic constants, piezoelectric coefficients, and di- 
electric constants can be written as 
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The variational formulation using curvilinear coordinate can be derived using the 
above relationships, and can be expressed in compact notation as 

° = i: {- 1 ( s [u ’ ds+ i{ s ^ ^ p *c< + 



( 118 ) 
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Discrete-Layer Theories 

To analyze thick laminated piezoelectric shells, discrete-layer theories are used in cre- 
ating the numerical models. The governing equations in matrix form are available for 
both cylindrical and general shells. 

Background And Motivation 

Most theories developed for analyzing laminated piezoelectric composites include limi- 
tations which may not represent the true behavior for certain applications. For the analy- 
sis of thick laminates under small deflection, the use of Kirchhoff hypothesis as theoretical 
background causes an undesirable approximation of both interlaminar and intralaminar 
stress components. In addition to this limitation, many previous studies have made use 
of an equivalent force representation of induced strain actuation in the piezoelectric lam- 
inate. This kind of approach does not solve the coupled equation of piezoelectricity. To 
analyze thick laminated piezoelectric shells, cylindrical or general, the above limitations 
need considerable refinement. This is the direction of this research. 

In this research, discrete-layer theories developed by Pauley [Si], Reddy, Barbero, 
Teply [32, 43, 44], Robbins [53], and Heyliger and Saravanos [59] are used. The main 
reason for using these approximations is that the coupled relationship between elastic 
and electric variables can be exactly represented. Additionally, the limitations from the 
Kirchhoff hypothesis can be avoided. 

In discrete-layer theories, the kinematic assumptions as shown in Figure 5 defined for 
the laminated plate element are : 

1. Through the thickness of composite elements, arbitrary variations are allowed within 
the in-plane displacement components and the electrostatic potential (i.e. u,v, and 
<f > ), and 

2. There are three assumed forms of the transverse displacement component ,w, which 
are used to develop two separate theories. 
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• A constant transverse displacement through the thickness, 

• A distributed transverse displacement which is identical in form to the other 
variable through the thickness. 

• A distributed transverse displacement of which the form is not completely 
layer-wise like the other variable through the thickness as shown in Figure 6. 


Within the laminated cylindrical shell element, the first assumption yields reference- 
surface displacement components and electrostatic potential that are assumed to have 
arbitrary variations through the thickness (r-direction for cylindrical shell, (^-direction 
for general shell). In the second assumption, the transverse displacement component, 
w(r, 0, z ) or u/(<f, 77, <£), is with the direction normal to the shell surface (r- or £- direction). 

The theory developed with a constant transverse displacement (i.e. w(r,Q, z) = w(Q, z) 
or u>(f , 77, () = u;(£, 77)) is for simpler and more economical in computation. However, poor 
accuracy is given when calculating the interlaminar stress or analyzing thick laminates. 
As with the case of the constant w formulation for the plate, it is impossible to capture 
the through-thickness actuation strain using this form of approximation. With the second 
assumed form (i.e. w is function of r,0,z or function of £,r],Q, the approach is more com- 
plex and more expensive, yet much more accurate. Both approaches will be investigated 
in this research. With the use of discrete-layer theories, the variational formulation is 
further derived for the type of approximation used in this approach. 


Governing Equation in Matrix Form for Cylindrical Shell 


In this section, a re-arranged variational formulation of governing equation is presented 
where a configuration of simultaneous equations is obtained. Secondly, an approximation 
of displacement and potential variables is given. The development of the discrete-layer 
theory for piezoelectric laminated shell begins with this. Finally, a matrix form is ex- 
pressed for further computational models. 
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Here each subscript behind the right brace shows the connection between the variational 
form and the matrix form that is going to be derived. 

The displacement and potential variables are approximated using linear combinations 
of the form 
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If a distributed transverse displacement is assumed in the discrete-layer theory, one- 
dimensional Lagrangian interpolation polynomials ipj(r) can be used for the through thick- 
ness approximation in the above (4>j{r) = V , j( r ) = x Pj J { r ) = V’/( r ) = ? /’j( r ))- Thus, n — 1 
is defined as the number of subdivisions through the thickness. For better results, the 
number of subdivisions should be larger or equal to the number of layers in the laminate. 
The in-surface approximation for the cylindrical shell is assumed using two-dimensional 
functions (0,z). The related number m is the total number of functions for the in-surface 
approximation. By replacing the variables with the approximations, the governing equa- 
tions can be expressed in matrix form as 
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Both the [M] and [K] matrices are symmetric, so, [A' 12 ] = [A" 2 i] r , [A' 13 ] = [A 3 i] r , 
(/<■«] = [X 41 f, [Xa] = [Xa] T , [X 24 ] = (X 42 ] T , and [X*] = |A'„] T . The structure of 
the submatrices ([A’hJJA'u],...,) is explained in Figure 7. Here the boxes drawn in the 
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submatrix locates the members of submatrix using certain approximate functions 'f'a and 
^6. Those box areas are noted as [/Cn]n, [ifi 2 ]i 3 > etc. The expressions in these cases are 
given as follows. 
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The sub- matrices [A,j], [jB.j], .... are expressed as 
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(152) 
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Here N is the total number of layers (iV = n — 1). 

However, if a constant transverse displacement is assumed, then = i>j{r) = 

^(r) = ipj(r) but r ) = constant, say /c, and several of the above equations must be 

changed. 
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dOdz 
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A Discrete-layer Element for General Shell 

To model the laminated shell with arbitrary geometry, a curvilinear representation of 
a discrete- layer shell element is developed using the kinematic assumptions of discrete- 
layer theories. This is shown schematically in Figure 8. Within the shell, continuity of 
the in-surface displacement components at a number of locations through the thickness 
(i.e. usually at the joint between layers) is required as shown in Figure 9. For the 
two different forms of the transverse displacement (out-of-surface displacement), different 
requirements need to be satisfied. If a distributed out-of-surface displacement is assumed, 
the same continuity of this displacement component is required as that of the in-surface 
components. However, if a constant out-of-surface displacement is assumed, the continuity 
of this displacement component is necessary at only one location through the thickness, 
such as the bottom or top surface of the shell. 

The weak form of the governing equation using curvilinear coordinates was shown ear- 
lier and is used to allow a much easier integration for the shell with arbitrary geometries. 
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To impose continuity across the interface of adjoining elements witli no loss of generality, 
the resulting degrees of freedom are maintained at the local level throughout the analy- 
sis. This type of formulation can give a much more accurate representation of adjoining 
element sides. With the defined Gauss points, the numerical integration can be accom- 
plished by mapping the discrete-layer shell element into a parent element that is similar 
to the discrete-layer plate element described for the plate earlier in this report. Quadratic 
in-surface approximations will be used in the calculations to allow representation of the 
curved edges. This order of approximation can be easily modified if necessary. 

To define the geometry of a laminated shell element and to determine the Jacobian 
matrices [J] and [J*] .matrices, the geometry of shell element [69] is described in standard 
fashion as 


\ 

X 

n f Xi 

< y 

> = 7,C)< yi > 

, 2 . 

1 J 


Here (x,j/,z) represents the location of a generic point within the shell element, and 
(xi,yi,Zi) are the locations of element nodes. The function N, indicates the shape function 
corresponding to the node i. If the curvilinear coordinate systems used in all layers 
are same, then the geometry of this element can be defined using only the locations of 
the nodal points on any two layer surfaces, usually the top and bottom surfaces of the 
element. However, if more than one curvilinear coordinate systems is required, more 
nodes are needed for the approximation. The displacement components ( u,v,w ) and 
electric potential in local system can be approximated in this fashion by the following 
equations, 
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The member Dij of the [jD]i 2 xi 2 matrix can be calculated by 
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By replacing the variables with the approximations, the governing equations can be rewrit- 
ten in matrix form. Let t/>“(C) = ^j(C) = V’j (C) = V’j(C) but V’J'(C) > s not necessarily equal 
to t/)j(C). The expressions of the [M\ and [K] for the general shell are given as follows. 
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Where the members of [A,j], [B tJ ] , ... , are expressed as 
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c?r = £ / Dk.*7*T<H 
1=1 

(220) 


(221) 


(222) 

V/ C,4 'n »*?,,, 

Cij tih Dkm 5C dc 

(223) 


If a distributed transverse displacement is assumed in the approximation, then ?/>j = if;™. 
Therefore, [A] = [B] = [C],[A] = [i] T = [B\ = [B] T = [C] = [C T ], and [A] = [5] = [C\. 
However, if a constant transverse displacement is assumed (t/>j 1 '(C) = constant), then 
[B] = [ B ] = 0 and [C] = [C] = [C] = 0. The fundamental behavior of the shell is shown 
for several layers in Figure 11. 

Static solutions to problems involving this type of approximation are fairly straight- 
forward, However, there are several difficulties involved in dynamic analyses. These are 
described separately below. 


Dynamic Analyses 

The kinetic energy of the system is involved when developing the Hamilton’s princi- 
ple into the matrix formulations. Hence, piezoelectric vibration analysis is possible using 
solutions of equation 128. The primary problem in this research to lincl the free vibra- 
tion characteristics (i.e. natural frequencies and corresponding mode shapes) for given 
arbitrary shells and cylindrical shells, as well as the behavior of shells under a forced 
vibration. 


Free Vibrations 


Since no external force term is assumed in free vibration problems, equation 128 can 
be rewritten as [59] 
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(226) 


{A} = 


{«} 
M 
. M , 


If periodic motion is assumed, equation 128 becomes an eigenvalue problem, expressed as 

{[/?]-u-MJV]}{A} = {0} (227) 

Here |/?j is formed using static condensation, which is necessary to eliminate the potential 
variables. This matrix is given by 

'/?] = [l< AA ] - [l< A *\ [K**] (228) 

Here u is the natural frequency and the corresponding {A} presents the mode shape. 
Forced Vibrations 


It is possible to analyze the behavior of a laminated piezoelectric shell under a forced 
vibration, especially a harmonic excitation. According to the governing equation in matrix 
form, there are several ways to cause excitation in the system : (1) applying traction 
forces {/(£)}, (2) applying an electric field {<f>{t)}, (3) increasing / decreasing charges in 
the surface electrode, (4) mix of (1) and (2), and (5) mix of (1) and (3). 

If the system is excited by surface traction forces which are harmonic, and no electric 
field is applied, equation 224 can be rewritten as 
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( 229 ) 


The equation of motion of this system can 
the equation of general laminates. 


{unlmotcn} j 
be expressed to the form that is the same as 


[M\ {A} + [A' a4 ] (A) = {/(()} (230) 


After the displacements, {A}, are solved, j {A} can be used to predict the surface 
charge created by shape change of the structure. The displacements can be discovered 
by measuring the surface charges which are collected through a surface electrode to the 
outside detector. Furthermore, if a static electric field is applied, the equation of motion 
becomes 


[A/] {A} + [A aA ] (A) = {/(()} - [A A ‘] M 
The surface charges can be predicted by solving {A} + {o}. 


( 231 ) 
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If the system is excited by an applied electric field that is changed harmonically, and 

no traction force is applied, equation 224 can be rewritten as 
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The equation of motion of this system can be written as 

[M\ {A} + [l< AA ] {A} = - \l< A *} {^(i)} (233) 

The equation of surface charges becomes {Q} = [/<f* A ] { A} + { o(Z ) } . 

The third way to excite the system is harmonically increasing / decreasing the charges 
that store in the surface electrode. Equation 224 can be rewritten as 
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Above equation can be expressed as 


[M] {A} + [/?] {A} = - [l< A +] [. K **}’ {Q(t)} 
[/?] = [k aa ] - [k a *] [k^Y' [a* a ] 

The corresponding voltage can be found using 

W=[A'«]"{[Q(i)l-[A' M ] (A}} 
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If the system is excited by traction forces and an applied electric field, which are both 
changed harmonically but not necessarily in the same frequency, equation 224 can be 
rewritten as 
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The equation of motion becomes 


[Ml {A} + [*“] (A) = (/(()} - [/< a *] W (239) 

The surface charges can be predicted using [2r* A | {A} + [ A ** | { <p( 2 )). Since the motion 
of the intelligent structure can be sensed, it is possible to apply a required electric field to 
create or change the damping mechanism of this structure system. This effect is proposed 
to control (fully or partially) the vibration in this structure. 
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If the system is excited by traction forces which is harmonic and increasing / decreasing 
the surface charges in the electrode harmonically, equation 224 can be rewritten as 
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The equation of motion becomes 
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[M\ {A} + [/?] {A} = {/(*)} - [A'^] [K**]' 1 {Q{t)} (241) 

The corresponding voltage can be found using equation 224. These equal ions also conclude 
that if the current that flows into the surface electrode and the voltage in the electrode are 
measured, the motion of this structure could be predicted. Furthermore, by increasing / 
decreasing the charges, vibration of the shells could be controlled. 1 hi.-> type of behavior 
can be used in the active tip-clearance control of engine blades. 


Computational Models 

A prototype discrete-layer shell element program has been completed and is included 
as part of this report. 
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6. EXACT SOLUTIONS 

Exact solutions are developed for predicting the coupled electromechanical vibration 
characteristics of simply-supported laminated piezoelectric plates. The three-dimensional 
equations of stress and electric displacement equilibrium are solved using the assumptions 
of the linear theory of piezoelectricity. The through- thickness distributions for the dis- 
placements and electrostatic potential are functions of eight constants for each layer of 
the laminate. Enforcing the continuity and surface conditions results in a linear system 
of equations representing the behavior of the complete laminate. The determinant of this 
system must be zero at a resonant frequency. The natural frequencies are found numeri- 
cally by first incrementally stepping through the frequency spectrum and refining the final 
frequencies using bisection. Representative frequencies and mode shapes are presented 
for a variety of lamination schemes and aspect ratios. 

This chapter is included because it represents the results of a tangential study under- 
taken as part of this research that was effectly a matter of necessity. Few solutions exist 
for laminated piezoelectric plates, and this phase of the study was crucial in determining 
the effectiveness of the discrete-layer approximations described as the main part of this 
research. 

Introduction 

The behavior of linear elastic laminated plates composed of dissimilar orthotropic ma- 
terials have been studied for a number of geometrical configurations, lamination schemes, 
and boundary conditions. The exact solutions of the equations of motion for these solids 
have only been obtained for the limited case of simple support. Static solutions have been 
comprehensively studied by Pagano [10, 11]. For the dynamic case, the two-dimensional 
case of cylindrical bending has been considered by Jones [76, 77] for two-layer cross-ply 
and angle-ply laminates, with the exact natural frequencies and mode shapes being ob- 
tained for a number of aspect ratios. The three-dimensional laminated plate geometry has 
been studied by Srinivas and coworkers [7S, 79] using an exact solution for a rectangular 
plate with simple support. 

Studies of the linear vibrations of finite laminated piezoelectric plates are limited. The 
monograph of Tiersten [SO] provides a comprehensive study of the governing equations, 
fundamental behavior, and exact and approximate solution methodologies for single-ply 
piezoelectric plates. Studies involving piezoelectric laminates have for the most part been 
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confined either to infinite plates or approximate theories and solution techniques. These 
include the studies of Pauley [81], Ricketts [82, 83], Lee [84], and Yong and coworkers 
[85]. 

The objective of the present study is to develop exact solutions for the static response 
and natural frequencies of free vibration for simply-supported, laminated, rectangular 
plates composed in part of orthotropic piezoelectric layers. The sirnplv-supported plate is 
one of the few geometries for which the in-plane functions can be selected to exactly satisfy 
the governing equations and boundary conditions. The through-thickness distributions 
are evaluated and found to be a function of eight unknown constants that are frequency 
dependent. Imposition of the interface and surface conditions for the complete laminate 
results in a linear system of equations that must be iteratively solved for each frequency. 
This study builds on the work of Pagano [11], who developed exact solutions for static 
behavior of elastic laminates. 

The developed solution and results presented here should provide a good basis for 
comparison for approximate plate theories. The results could be especially important 
for thick plates, in which case many theories provide poor approximations, and will help 
establish the limitations and ranges of applicability of other approximate methods. 

Exact solution 

Governing Equations 

The geometrical configuration of the laminate is such that that the thickness dimension 
of the laminate coincides with the z-direction, with the lengths of the plate in the x and y 
directions denoted as L r and L y , respectively. Each layer of the laminate can have elastic, 
piezoelectric, or conducting material properties. The general problem considered in this 
study is to determine the behavior of the elastic and electric field components throughout 
the laminate under periodic vibration with arbitrary surface conditions. 

The constitutive equations for each layer in the laminate are assumed to be those of 
a piezoelectric material. These equations can also be used to represent purely elastic or 
conducting layers, and are given by [80] 

< 7 , = Co-5, - tfiEj (242) 

Di = eijSj + tijEj 
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Here <j x - are the components of the stress tensor, C x; - are the elastic stiffness components, 
S j are the components of infinitesimal strain, e JX - are the piezoelectric coefficients, E x are 
the components of the electric field, D; are the components of the electric displacement, 
and are the dielectric constants. The standard contracted notation has been used here. 
In this study, the poling direction is coincident with the x 3 or z axis. 

The displacement components u x , where Ui = u, U 2 = v, and u 3 = w, are related to 
the strain components through the relations 


c - ^ i\ 

ij 2 [dxj dxi) 

To be consistent with Equation 242, the conventional notation for the strain indices has 
been used, i.e. Sn = Si, S 23 = S^, etc. The electric field conipoiienis can be related to 
the electrostatic potential (f> using the relation 




(244) 


For the materials used in this study, it is assumed that the non-zero components of the 
rotated piezoelectric tensor e,j are e 3 i, e 32 , 633 , e 2 ^, and eis. The elastic stiffnesses C tJ are 
those of an orthotropic material, and the dielectric constants are given by tu, 022 , a.nd 


£33- 

The stress equations of motion are given by 


a ij,j = pui 

and the charge equation of electrostatics is given as 
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D u = 0 (246) 

Substituting the constitutive relations, the stress-strain relations, and the field-potential 
relations into Eqs. 245 and 246 gives the governing equations of the laminae in terms of 
the displacement components u, v and w and the electrostatic potential <j> as 
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These represent the coupled governing equations for a single layer within the laminate. 

For the problems considered in this study, an arbitrary number of laminae are assumed 
to be perfectly bonded together. At the top and bottom surfaces of the laminate, a 
given load, displacement, potential, or charge can be specified. A number of surface 
conditions could be treated using this methodology, but the cases of primary interest are 
those in which both upper and lower surfaces are traction- free. It is also of interest to 
consider the influence of electric surface conditions, and the cases of specified homogeneous 
potential and transverse electric displacement are treated here. The laminate is assumed 
to be simply supported, and the vertical edges of the laminate are assumed to be fixed 
at zero (grounded) potential. Hence along a plate edge, the normal stress, tangential 
displacement, transverse displacement, and electrostatic potential are specified to be zero 
regardless of the remaining conditions on the remaining laminate surfaces. 

At each interface between layers, continuity conditions of displacement, traction, po- 
tential, and electric displacement must be enforced. Using an indexing scheme, the con- 
ditions for the i-th layer can be expressed as, for example, 
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Here i represents the layer number, with i=l the top layer, each layer has an individual 
coordinate system with the origin at the left end in the center of the layer, and h, is the 
thickness of the i-th layer. Similar interface conditions exist for v, w, o. a z , t X j, r yz , and 
D z . At a single interface of a laminate with n plies, there are six conditions related to the 
elastic variables and two conditions related to the electrostatic variables for a total of 8(n- 
1) interface conditions for the complete laminate. At both the top and bottom surfaces, 
there are three elastic boundary conditions and one electric condition for a total of 8 
conditions. The surface conditions require specification of one variable from each of the 
pairs (u,t X2 ), (v,r y2 ), (w,cr 2 ), and (^>,D 2 ). Enforcing all conditions leads to 8n equations 
relating the variables within all layers of the laminate. 

Method of Solution 

Solutions for the displacement components and the electrostatic potential are sought 
in the form 


u(x,y,z,t ) = U(z) exp(ju>t) cos pi sin qy = U exp(jut) exp(sz) cos par sin qy (252) 
v(x,y,z,t) = V(z,exp(jut)t) sinpx cos qy = V exp(jut) exp(sr) sin px cos qy 
w(x,y,z,t ) = W (z , t) exp(j ut) sin px sin qy = W exp(jut) exp(sz) sin pa- sin qy 
<f>(x,y,z,t ) = $( 2 , t)exp(jut) sin px sin qy = $ exp( jut) exp(sz) sin pxs'm qy 

Here the overbarred terms are constants, s is an unknown number, p = m tt/L*, q = n 
7r/L y , and z is the local layer coordinate whose origin is at the center of each laminae. 
Substitution of these expressions into the equations of motion and the charge equation 
yields the system of equations 
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The Aij elements of this matrix are 
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/1 22 _ Ceep 2 + C 229 2 — P u 1 A 23 — —q(C23 + C 44 ) A 24 — — <jf( e 32 + e 24) (254) 


A33 — C55P 2 + C445 2 — pu; 2 A 34 — eis p 2 + e 24<? 2 A44 = ~£n/' 2 — ^22^ 2 

A non-trivial solution to this system requires a zero determinant. The resulting char- 
acteristic equation is 


As 8 + Bs 6 + Cs 4 + Ds 2 + E = 0 (255) 

Expressions for the coefficients of this polynomial are given in the Appendix. This equa- 
tion can oe written as the fourth-order equation 


r 4 + cr 3 -f- dr 2 + tr + / = 0 

where 


(256) 
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The roots of this equation a.ve a function of the material properties, the laminae geometry, 
and the frequency u>. The roots can be real, imaginary, or complex, which results in 
different forms for the solutions for each variable. Regardless of the nature of the roots, 
the solutions for a given value of s are based on the original form for the solution for U(z). 
The remaining components can then be computed using Eq. 253, which is rearranged in 
terms of the unknown constants as 
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General expressions that can be used to evaluate the constants V, IT, arid <j) are con- 
structed as a function of the real, imaginary, or complex roots. These are 


V(s) = lllfl+JlllLlIllO 
K ’ D(s) 


( 260 ) 
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( 261 ) 


Hs) = W) 

D{s) — d\S 6 C^S ^ “t" C?3.S 2 “t“ C?4 

The constants d,- and f,j are defined in the Appendix, and will change in nature depending 
on the type of root. These values and the corresponding solutions for the elastic and 
electric field components corresponding to each type of root are developed separately 
below. 

Case 1: Real roots for r 

Given n real roots for r, the 2n roots for s can be obtained using equal ion 256. These 
roots are either real or imaginary depending on the sign of r. Following the nomenclature 
used in Pagano [10] and Heyliger and Brooks [S7], the solution for the displacement 
components and electrostatic potential corresponding to the these roots can be written 
in either case as 


(262) 

(263) 


U(*) = 't Ui (z) V(z) = ■£ LjUj(z) W(z) = f: M,w,{z) 

j=l j = 1 j= 1 

where 


<!>(:) = 

j= 1 


Ui = FjCj{z) + GiSj(z) (265) 

Wj = GjCjlz) + ajFjSjiz) 

Here Fj and Gj are real constants, there is no summation on j, and the functions C and 
S and the values mj and otj are defined as 


Cj(z) = cosh(mjz) 

Sj(z) = sinh(mjz) a, =1 (r X 0) 

(266) 

Cj(z ) = cos (rrijz) 

Sj(z) = sin(mjz) aj = — 1 (»’ -< 0) 

(267) 


m, = |sj| 

(268) 
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The coefficients L j, Mj, and Nj are more specific representations of the parameters given 
in Eqs. 269-271, and are given in this case as 


Lj = (/nmj + fuoijm] + / 13 ) (269) 

Mj = (fnrrij + hiOtjTn) + /») (270) 

iVf = ~ + hiotj-m) + fjz) (271) 

where the determinant Dj is given by 

Dj = diQj-m® + d 2 mj + d 3 a_,mj + d A (272) 


Using the constitutive equations in (242), the corresponding expressions for the stress 
and electric displacement can be computed as 


<7, = s'm px sin qy^2[— pC, i — yC.iLj 4- (273) 
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Di = sinpxsin<?y£[-e3ip - tyiqLj + (277) 

i= i 

771 ^ TTl^ 

+ MrnjCtj + fa 3 ) + (330 31m* + f^rn^aj + J:n)]Wj(z) 

Here i=l,2,3 corresponds to x,y, and z for the stress and electric displacement components. 

A special case in which there are real roots for r is the non-piezoelectric elastic layer. 
The e,j=0 in this case and the elastic and electric fields uncouple. The elastic solution 
has been given by Pagano [10], and the results are not repeated here except to note that 
the elastic field behavior is represented by six roots and six unknown constants within 
the layer. This corresponds to the 6 interface/boundary conditions (three displacements 
and the 0 { z stress components) for a single layer. The electrostatic behavior in this case 
is represented using the two roots 
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The potential and transverse electric displacement components in this case are given by 


2 

<j>(x,y,z) = sinpxsmqy Bjexp(rijz) (279) 

j=i 

2 

D z = —€33 sin px sin qy ^ Bjrij exp(njz) (280) 

i= 1 

Case 2: Complex roots for r 

The elastic, electric, and geometric properties for some laminae yield complex roots. 
These appear in conjugate pairs, which result in the final roots for s in the form ±(a ± ib), 
where i = y/—l and a and b are positive constants. The solution for U(z) corresponding 
to these roots can be expressed as 


U(z) = Cie az cos bz + C 2 e az sin bz + c^e az cos bz + c*e sin 6 -r (281) 

where ci,...,c 4 are real constants. Following some algebraic manipulations and using Eqs. 
(260-262), the solution for V(z) can be expressed as 
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V(z) = c 1 e°*(ri cos bz — fi; sin bz ) + C 2 e“*(n 1 cos bz + T\ sin 6~) + 
c 3 e -02 (ri cos bz -f sin bz) + c 4 e~ az (~Cli cos bz + T i sin bz) 


(282) 


Here Ti = + i6)] and f2i = S[V(a + i6)]. Similarly, the final expression for W(z) 

can be expressed as 


W(z) = cie az [(aT 2 — 6Q 2 )cos6z + (— 6T 2 — afi 2 )sin6z] + (283) 

c 2 e az [(6T 2 + afi 2 ) cos bz + (aT 2 — 6fi 2 ) sin bz] + 
c 3 e~ az [(6fi 2 — aT 2 ) cos bz + (— 617 2 — afi 2 ) sin bz] + 
c 4 e~ az [(6r 2 -f aH 2 ) cos bz + (— aT 2 + 6Q 2 ) sin 6c] 

where T 2 = 3£[iy(a + i6)] and fi 2 = £r[M (a + i6)]. The final expression lor o can be 
obtained in similar fashion to yield 


( p(z ) = C\e az [(aT 3 — bCl 3 )cosbz + (— 6T 3 — afl 3 )sin6z] + (284) 

c 2 e az [(6T 3 + aff 3 ) cos bz 4- (aT 3 — 6Q 3 ) sin 6c] + 
c 3 e -az [(6fi 3 — aT 3 ) cos bz + (— 6T 3 — <zfi 3 ) sin 6c] + 
c 4 e -oz [(6r 3 + afl 3 )cos6z + (— nT 3 + 6n 3 )sin6c] 

where T 3 = 3£[</>(a + i6)] and T 3 = £r[^>(a + t"6)]. 

The functions for displacement and potential in each case must be combined with 
the solutions corresponding to the remaining roots to construct the complete solution 
for a given layer. The expressions for the stress and electric displacement components 
can be obtained by the appropriate differentiation and combination with the constitutive 
equations as given in Eq. 242. Because of their length and relative ease of calculation, 
these are not given here. 

Solution for the Laminate 

The elastic and electric field components within each layer are expressed in terms of 8 
unknown constants. These are determined using the interface and continuity conditions 
at the upper and lower surfaces of each lamina. For an elastic/dieloclric layer, there 
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are 6 constants corresponding to the elastic part of the solution and 2 constants for the 
electrostatic solution. For most common piezoelectric materials, however, the roots fall 
into one of two categories. For some piezoelectric materials, the four roots for r are real. 
Correspondingly, there are 8 constants (Fi,...,F 4 and Gi,...,G 4 in Eq. 265) that uniquely 
define the fields within the layer. For other piezoelectric materials, there are 2 real or 
imaginary roots and 2 complex conjugate roots for r. Hence the eight constants are F t , 
F 2 , Gi, and G 2 from Eq. 265 and Ci,...,c 4 from Eqs. (281-284). Following the solution 
of the total system of equations for the constants, the solution for any component can be 
computed at any location within the laminate. 

For static solutions, the coefficient matrices are all known as are the right-hand side 
elements. The solution for the unknown coefficients reduces to solving a linear system of 
equations. 

A resonant frequency u is unknwon a priori , yet the roots and through-thickness field 
distributions are a function of this value. An iterative scheme was used 10 evaluate these 
parameters. For the vibration problem, a necessary and sufficient condition for a non- 
trivial solution to exist is that the determinant /? of the coefficient matrix [A] multiplying 
the unknown constants is zero. The [A] matrix contains the final coefficients multiplying 
the unknown constants following imposition of all interface and surface conditions and is 
written as 


[A]{A} = {0} (285) 

where A is the column vector of coefficients. The zero determinant requirement establishes 
the necessary relationships for the calculation of the natural frequencies and the modal 
shapes for the plate. 

To evaluate the roots of the characteristic equation, the following procedure was used. 
First, the frequency was stepped through a sequence of frequencies that are an increment 
of the lowest expected frequency. This can either be estimated or approximated from 
other plate theories. The sign of the determinant is then computed for each value aiid 
recorded. Once a sufficient number of sign changes have been noted, bisection was used 
to refine the values of the true frequencies using the sign-change values as the bounding 
initial guesses. 

It can be a difficult numerical problem to evaluate the determinant of the matrix as 
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it becomes singular, especially if computed on a limited-precision machine. In this study, 
the determinant 0 was computed using 


0 = IIJLjAi (286) 

where A,- is the i-th eigenvalue of the coefficient matrix [A]. This matrix is square and of 
order m, where m=8n and n is the number of layers within the laminate. The eigenvalues 
were computed using the QR algorithm in extended precision. Attempts to compute the 
eigenvalues using the diagonals from an LDU decomposition and QR with di able precision 
failed. 

Following computation of the frequencies uj , the constants corresponding to these 
frequencies that define the modal displacements and potential are readily computed by 
solving the system of equations with the first constant in A arbitrarily set to one. The 
reduced [A] matrix is no longer singular after this step, and the relative values for the 
remaining constants can be found. The through-thickness mode shapes can then be 
calculated using these constants. 

Clearly, a different characteristic equation can be generated in determinant form for 
each combination of p and q in Eq. 259. These describe the in-plane modal characteristics 
of the elastic and electric fields. For each of these equations, there are an infinite num- 
ber of eigenvalues corresponding to the through-thickness modes. The cases of interest 
in this study were the thickness modes corresponding to the fundamental in-plane mode 
with m=n=l in Eq. 252. Of particular interest is the relationship between the displace- 
ment components and the electrostatic potential, as this information could be used to 
electrically excite or sense the various modes of vibration. 
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7. NUMERICAL EXAMPLES AND RESULTS 


This section contains results of applications of the computational and analytic models 
developed as part of this study. This includes 

• Results from laminated piezoelectric plate element. 

• Semi-analytic solutions for plates. 

• Results of exact solutions for simply-supported plates. 

The shell elements were in the process of being tested at the time of writing this report. 
Several simple test cases were being studied (simple electric field, uniform strain applica- 
tion, etc.). One of the difficulties of benchmarking these elements is the lack of solutions 
available for laminated piezoelectric shells. Although a proposed area for future study is 
development of an exact solution for this class of problem, no such solutions are currently 
available. Because of the preliminary status of the results for the shell, they were not 
included in this report. 

The Plate Element 

In this section, several examples are investigated in order to test the accuracy, validity, 
and range of applicability of the theories presented in section 4. Initial examples are for 
the purely elastic plate, with later examples to include piezoelectric layers. 

Simply Supported Cross-Ply Plate Under Sinusoidal loading 

A simply supported symmetric cross-ply (0/90/0) laminated rectangular plate under 
a sinusoidal transverse distributed load on the upper surface is considered. This example 
is chosen because it has an elasticity solution developed by Pagano [10] and therefore, 
provides an excellent tool to check the level of accuracy of the two theories presented 
herein for purely elastic laminates. 

Each of the three layers has equal thickness and is idealized as a homogeneous or- 
thotropic material with the following properties: 

E x = 25* 10 6 psi 
G 23 = 0.2*10 6 p,sz 
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E? = £*3 = 1.0 * lO^psi 
G \2 = G13 = 0.5 * 10 


V\2 — ^23 — ^13 “ 0.25 


The domain of the rectangular plate is0<x<a,0<j/<6, and 0 < z < H with 
b = 3a. The transverse sinusoidal distributed load is given by: 

7TX TTV 

q(x,y) = q 0 sin( — )«'n(-) (287) 

a b 

Because of the symmetry of the problem, only a quarter of the plate is modeled 
(0 < x < |, 0 < x < |) and the corresponding boundary conditions are given by: 


u(x, 0, z) = 0 

u ( a / 2 iy, z ) = 0 
u(0,y,z) = 0 
v{x,b/2,z) = 0 

it’(0,y,z) = 0 

u;(x,0,z) = 0 


(2SS) 


The deflections and stresses are computed by using two theories, o nr considering w 
varying linearly through the thickness of the plate and the other one. with w constant 
through the thickness. The stresses are determined using the computed displacements 
and the constitutive relations 38 and 39, and they are evaluated at the closest gauss 
point to the the following locations: 

a t (a/2,b/2,H) 

0y{a/2,b/2,H/Z) 

rxy(0,0,/O 

Tx*(0, 6/2,0) 

r yz (a/ 2,0,0) (2S9) 

Regarding the transverse displacement iu, it is evaluated at the center of the rectan- 
gular plate. 

The fortran programs written to compute the deflections and stresses according to the 
two theories presented herein were run for four different meshes : a CxG linear element mesh 
with 6 and 12 layers in the z-direction and denoted by L6x6-6 and L6x6-12 respectively, 
and a 3x3 quadratic element mesh with 6 and 12 layers in the z-direction and denoted 
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Mesh 

a/h=2 

a/h =4 

a/h=10 

a/h=20 

a/h=5U 

a/h= 100 


10"° 

10"° 

lO-* 

10"* 

10 -J 

10-*' 

L6x6-12 (F) 

0.64710 

0.17956 

0.90463 

0.46290 

0.49566 

0.22721 

L6x6-12 (FR) 

0.65119 

0.18103 

0.91918 

0.48600 

0.64715 

0.50479 

L6x6-12 (R) 

0.65163 

0.18128 

0.92210 

0.48831 

0.65075 

0.50767 

L6x6-6 (F) 

0.62983 

0.17752 

0.90058 

0.46212 

0.49571 

0.22719 

L6x6-6 (FR) 

0.6 

0.17894 

0.91494 

0.48513 

0.64688 

0.50470 

L6x6-6 (R) 

0.63405 

0.179192 

0.91786 

0.48743 

0.65048 

0.50759 

Q3x3-12 (F) 

0.64553 

0.17968 

0.91673 

0.48660 

0.64773 

0.50256 

Q3x3-12 (R) 

0.64570 

0.179699 

0.91695 

0.48711 

0.65024 

U. 50742 

Q3x3-6 (F) 

0.62823 

0.17762 

0.91255 

0.48573 

0.64744 

0.50245 

Q3x3-6 (R) 

0.62841 

0.17765 

0.91278 

0.48625 

0.64998 

0.50733 

EXACT 

0.65327 

0.18055 

0.91891 

0.48763 

0.65060 

0.50766 


Table 1: Transverse Displacements at the Center of Plate, w Varying Linearly Through 
the Thickness 


by Q3x3-6 and Q3x3-12 respectively. The results are obtained using the full integration 
(denoted by F), the reduced integration (denoted by R), and the full integration for the 
bending terms and the reduced integration for the shear terms (denoted by FR). Tables 1 
and 2 shows the non-dimensionalized transverse displacements at the center of the plate 
for both theories and for a wide range of thickness/lenght ratio. 

The results for the case where a/h=4 are represented in figures 19 through 26 

Single Layer of PVDF Under Sinusoidal Load and Sinusoidal Polenlial 

A single square layer of PVDF is considered under two types of loads in using the 
two theories used in the previous section. First, a transverse sinusoidal load, q{x,y) = 
q 0 sin(^)sin(^) where q 0 = 1, is applied on the upper surface. Because of the symmetry 
of the problem, only one quarter or the plate is modelled and the boudary conditions used 
are as follows (see figure 17): 


u(x, 0, z) = 0 
t>(0,y,z) = 0 
™(0,y,*) = 0 
<p(x,0,z) = 0 
V>(x,y,0) = 0 


u(a/2,y,z) = 0 
v(x, 6/2, z) = 0 
w(x, 0, z) = 0 
<^(0,t/,2) = 0 
y(x,y,h) = 0 


The results for the case where a/h = 10 are represented in figures 27 through 37. 
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Mesh 

a/h=2 

a/h=4 

a/h=10 

a/h=20 

a/h=50 

a/h= 100 


io-° 

10-° 

: fir 3 

“TcF 5 

10" J 

lo--' 

L6x6-12 (F) 

0.66942 

0.18084 

0.90358 

0.46131 

0.49409 

0.22671 

L6x6-12 (FR) 

0.67351 

0.18232 

0.91807 

0.48423 

0.64409 

0.50232 

L6x6-12 (R) 

0.67421 

0.18226 

0.92098 

0.48652 

0.64766 

0.50518 

L6x6-6 (F) 

0.65305 

0.17883 

0.89952 

0.46055 

0.49397 

0.22670 

L6x6-6 (FR) 

0.65682 

0.18026 

0.91383 

0.48338 

0.64388 

0.50228 

L6x6-6 (R) 

0.65755 

0.18053 

0.91674 

0.48566 

0.64745 

0.50513 

Q3x3-12 (F) 

0.66766 

0.18096 

0.91557 

0.48480 

0.64466 

0.50013 

Q3x3-12 (R) 

0.66714 

0.18087 

0.91571 

0.48539 

0.64716 

0.50-192 

Q3x3-6 (F) 

0.65126 

0.17893 

0.91139 

0.48396 

0.64443 

0.50004 

Q3x3-6 (R) 

0.65072 

0.17884 

0.91153 

0.48447 

0.64695 

0.50488 

EXACT 

0.65327 

0.18055 

0.91891 

0.48763 

0.65060 

0.5766 


Table 2: Transverse Displacements at the Center of Plate, w . Constant Through the 
Thickness 


In the second case when the sinusoidal potential, y) = q 0 sin( ^)sin(^) where 
c p 0 = 1, is applied to the single layer of PVDF, the boundary condition are as follows: 


u(x, 0, z) = 0 
u(0 ,y,z) = 0 
™(0 ,y,z) = 0 
(p(x, 0, z) = 0 
V?(x,t/,0) = 0 


u(a/2,y,z) = 0 
v(x, 6/2, z) = 0 
u>(x, 0, z) = 0 
<p(0,2/,*) = 0 


The results for this case with aspect ratio equal to 10 ,a/h = 10. arc presented in in 
figures 3S through 4S. 

Single Layer of PZT4 Under Sinusoidal Load and Sinusoidal Potential 

A single layer of PZT4 is considered under the same types of loads as those examined 
with the single layer of PVDF. The geometry is shown in figure 17 wich is the same used 
for the example on section .The boundary conditions for the two cases, sinusoidal load 
and sinusoidal potential ,are the same as those in section respectively. 

The results for both cases, sinusoidal load load and sinusoidal poleulial, are shown in 
figures 49 through 56 and in figures 57 through 64. 

Four Layer PVDF/Graphite/Epoxy/PVDF laminate 

A simple supported square laminate composed of 4 layers, is considered in this section. 


Mesh 


< 

jy 

Try 

r it 


L6x6- 12(F) 

8.0302 

-6.0920 

0.8767 

-1.0261 

-0.2166 

0.2112 

0.5160 

0.1229 

L6x6-12(FR) 

8.0747 

-6.1332 

0.8849 

-0.10341 

-0.2183 

0.2129 

0.5066 

0.1233 

L6x6-12(R) 

8.1171 

-6.1494 

0.8866 

-1.0373 

-0.2191 

0.2134 

0.5072 

0.1230 

L6x6-6(F) 

7.1797 

-5.4417 

0.8313 

-0.9899 

-0.2054 

0.2021 

0.4940 

0.1063 

L6x6-6(FR) 


. 

. 






L6x6-6(R) 

7.2519 

-5.4885 

0.8403 

-1.0003 

-0.2075 

0.2041 

0.4979 

0.1061 

Q3x3-12(F) 

8.0676 

-6.1373 

0.8801 

-1.0324 

-0.2179 

0.2126 

0.-5049 

0.1238 

Q3x3-12(R) 

8.0693 

-6.1397 

0.8811 

-1.0334 

-0.2180 

0.2126 

0.5050 

0.1238 

Q3x3-6(F) 

7.2169 

-5.4850 

0.8350 

-0.9964 

-0.2065 

0.203-1 

0,1956 

0.1072 

Q3x3-6(R) 

7.2180 

-548645 

0.8354 

-0.9970 

-0.2066 

0.2034 

0.4957 

0.1072 

EXACT 

8.5323 

-6.4962 

0.9181 

-1.0709 

-0.2255 

0.2193 

0.5182 

0.1272 

L6x6-12(F) 

17.7195 

-17.0093 

1.6819 

-1.8451 

-0.4182 

0.4371 

1.3785 

0.1224 

L6x6-12(FR) 

17.8784 

-17.1650 

1.7038 

-1.8670 

-0.4233 

0.4421 

1.3927 

0.1199 

L6x6-12(R) 

17.9508 

-17.2284 

1.7070 

-1.8719 

-0.4242 

0.4432 

1.3941 

0.1196 

L6x6-6(R) 

17.1557 

-16.4718 

1.6430 

-1.8087 

-0.4107 

0.4298 

1.3756 

0.1088 

L6x6-6(FR) 

17.3022 

-16.6154 

1.6642 

-1.8299 

-0.4155 

0.4347 

1.3895 

0.1060 

L6x6-6(R) 

17.3748 

-16.6793 

1.6672 

-1.8346 

-0.4165 

0.4358 

1.3910 

0.1058 

Q3x3-12(F) 

17.8687 

-17.1653 

1.6911 

-1.8574 

-0.4218 

0.4407 

1.3865 

0.1263 

Q3x3-12(R) 

17.8790 

-17.1780 

1.6948 

-1.8611 

-0.4220 

0.4407 

1.3867 

0.1263 

Q3x3-6(F) 


-16.6279 

1.6529 

-1.8218 

-0.4141 

0.4333 

1.3833 

0.1126 

Q3x3-6(R) 


-16.6356 

1.6559 

-1.8244 

-0.4143 

0.4334 

1.3835 

0.1126 

EXACT 

18.3088 

-17.5873 

1.7405 

-1.9086 

-0.4301 

0.4490 

1.4076 

0.1291 

L6x6- 12(F) 

70.0011 

-69.9408 

4.0293 

-4.1974 

-1.1556 

1.1846 

4.0714 

0.1246 

L6x6-12(FR) 

71.2801 

-71.2188 

4.1162 

-4.2843 

-1.1794 

1.2084 

4.1527 

0.1014 

L6x6-12(R) 

71.6333 

-71.5713 

4.1301 

-4.3000 

-1.1839 

1.2133 

4.1535 

0.1013 

L6x6-6(F) 

69.5283 

-69.4711 

4.0047 

-4.1732 

-1.1491 

1.1781 

4.0688 

0.1126 

L6x6-6(FR) 

70.7863 

-70.7283 

4.0906- 

-4.2591 

-1.1725 

1.2016 

4.1498 

0.0892 

L6x6-6(R) 

71.1398 

-71.0810 

4.1044 

-4.2748 

-1.1771 

1.2064 

4.1506 

0.0891 

Q3x3-12(F) 

71.4377 

-71.3816 

4.0872 

-4.2585 

-1.1802 

1.2093 

4.1476 

0.1466 

Q3x3-12(R) 

71.5426 

-71.4925 

4.1117 

-4.2821 

-1.1806 

1.2095 

4.1486 

0.1471 

Q3x3-6(F) 

70.9790 

-70.9259 

4.0635 

-4.2351 

-1.1734 

1.2025 

4.1443 

0.1341 

Q3x3-6(R) 

71.0574 

-71.0074 

4.0880 

-4.2585 

-1.1737 

1.2027 

4.1454 

0.1347 

EXACT 

72.5977 

-72.5396 

4.1757 

-4.3490 

-1.1977 

1.2271 

4.2018 

0.1490 


Table 3: Stresses on Cross-Ply. w varying through thickness 
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a/h 

Mesh 

& X 

a y 


frz 


20 

L6x6- 12(F) 

0.2421 

-0.2422 

0.1102 

-0.1119 

-0.3442 

0.3472 

0.8049 

0.1187 

L6x6-12(FR) 

0.2553 

-0.2553 

0.1164 

-0.1181 

-0.3632 

0.3662 

0.8407 

0.0165 

L6x6-12(R) 

0.2567 

-0.2567 

0.1169 

-0.1186 

-0.3651 

0.3682 

0.8406 

0.0173 

L6x6-6(F) 

0.2418 

-0.2418 

0.1105 

-0.1122 

-0.3435 

0.3466 ; 

0.8040 

0.1029 

L6x6-6(FR) 

0.2548 

-0.2549 

0.1167 

-0.1184 

-0.3625 

0.3656 | 

0.839 i 

0.0332 

L6x6-6(R) 

0.2563 

-0.2563 

0.1173 

-0.1190 

-0.3644 

0.3675 | 

0.8396 

0.0340 

Q3x3- 12(F) 

0.2559 

-0.2560 

0.1155 

-0.1172 

-0.3646 

0.3676 1 

0.8575 

0.2314 

Q3x3-12(R) 

0.2567 

-0.2567 

0.1167 

-0.1184 

-0.3648 

0.367s : 

U.8577 

0.2331 

Q3x3-6(F) 

0.2556 

-0.2556 

0.1157 

-0.1175 

-0.3639 

0.3669 0.8565 

0.2135 

Q3x3-6(R) 

0.2562 

-0.2563 

0.1171 

-0.1188 

-0.3641 

0.3671 U.6566 ] 

0.2158 

EXACT 

0.2600 

-0.2601 

0.1178 

-0.1195 

-0.3696 

0.3727 . 0.8686 

0.2344 

50 

l 

L6x6- 12(F) 

1.1741 

-1.1740 

0.4864 

-0.4881 

-0.1577 

0.1580 j 1.5278 


L6x6-12(FR) 

1.5405 

-1.5406 

0.6387 

-0.6404 

-0.2070 

0.2073 ; 1.7991 


L6x6-12(R) 

1.5493 

-1.5493 

0.6421 

-0.6439 

-0.2082 

0.2085 i 1.7969 


L6x6-6(F) 

1.1740 

-1.1741 

0.4998 

-0.4915 

-0.1577 

0.1580 j 1.5257 


L6x6-6(FR) 

1.5440 

-1.5405 

0.6431 

-0.6448 

-0.2069 

0.2073 ; 1.7962 

. 

L6x6-6(R) 

1.5492 

-1.5493 

0.6466 

-0.6483 

-0.2081 

0.2084 | 1.7939 


Q3x3- 12(F) 

1.5373 

-1.5374 

0.6314 

-0.6331 

-0.2076 

0.2079 ! 2.1727 

0.5277 

Q3x3-12(R) 

1.5497 

-1.5498 

0.6417 

-0.6434 

-0.2082 

0.2085 i 2.1654 

0.5395 

Q3x3-6(F) 

1.5374 

-1.5375 

0.6358 

-0.6375 

-0.2075 

0.2078 ! 2.1696 

0.4888 

Q3x3-6(R) 

1.5496 

-1.5497 

0.6461 

-0.6479 

-0.2081 

0.2085 2 1624 

0.5014 

EXACT 

1.5690 

-1.5691 

0.6457 

-0.6475 

-0.2109 

0.2112 2.1927 

0.5393 

i 

100 

L6x6- 12(F) 

2.7435 

-2.7436 

1.1199 

-1.1216 

-0.3653 

0.3657 ! 1.2*175 


L6x6-12(FR) 

6.1300 

-6.1301 

2.5033 

-2.5050 

-0.8166 

0.8168 

; 1.2881 


L6x6-12(R) 

6.1651 

-6.1652 

2.5173 

-2.5100 

-0.8212 

0.8215 

1 .2702 


L6x6-6(F) 

2.7446 

-2.7446 

1.1285 

-1.1302 

-0.3653 

0.3656 

1.3565 


L6x6-6(FR) 

6.1312 

-6.1313 

2.5223 

-2.5240 

-0.8164 

0.8167 

1 .2825 


L6x6-6(R) 

6.1664 

-6.1665 

2.5365 

-2.5382 

-0.8211 

0.8214 

1.2646 


Q3x3- 12(F) 

6.0621 

-6.0622 

2.4590 

-2.4607 

-0.8146 

0.8149 

4.3995 

1.0025 

Q3x3-12(R) 

6.1668 

-6.1669 

2.5175 

-2.5192 

-0.8214 

0.8217 

4.3383 

1.0709 

Q3x3-6(F) 

6.0634 

-6.0635 

2.4778 

-2.4795 

-0.8144 

0.8147 

4.3936 

0.9265 

Q3x3-6(R) 

6.1681 

-6.1682 

2.5366 

-2.5383 

-0.8213 

0.8216 

! 1.3320 

0.9975 

EXACT 

6.2435 

-6.2436 

2.5308 

-2.5326 

-0.8319 

0.8322 

! 1.3916 

1.0652 


Table 4: Stresses on Cross-Ply. w varying through thickness 
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a/h 

Mesh 

&X 

d y 

Cry 

1 

V 


L6x6-12(F) 

6.8657 

-6.8657 

0.9860 

-0.9860 

-0.2220 

0.2220 

0.5082 

0.1205 


L6x6-12(FR) 

6.9080 

-6.9080 

0.9945 

-0.9945 

-0.2238 

0.2238 

0.5115 

0.1209 


L6x6-12(R) 

6.9345 

-6.9345 

0.9972 

-0.9972 

-0.2244 

0.2244 

0.5122 

0.1205 


L6x6-6(F) 

6.1660 

-6.1660 

0.9588 

-0.9588 

-0.2120 

0.2120 

0.4957 

0.1071 

2 

L6x6-6(FR) 

6.1960 

-6.1960 

0.9667 

-0.9667 

-0.2136 

0.2136 

0.4987 

0.1072 


L6x6-6(R) 

6.2231 

-6.2231 

0.9694 

-0.9694 

-0.2142 

0.2142 

0.4995 

0.1069 


Q3x3-12(F) 

6.9088 

-6.9088 

0.9909 

-0.9909 

-0.2233 

0.2233 

0.5099 

0.1215 


Q3x3-12(R) 

6.9116 

-6.9116 

0.9917 

-0.9917 

-0.2233 

0.2233 

0.5099 

0.1217 


Q3x3-6(F) 

6.2077 

-6.2077 

0.9638 

-0.9638 

-0.2132 

0.2132 

0.4974 

0.1081 


Q3x3-6(R) 

6.2092 

-6.2092 

0.9643 

-0.9643 

-0.2132 

0.2132 

0.4974 

0.1083 


EXACT 

1.62 

2.13 

0.268 

0.230 

0.0548 

0.0561 

0.257 

0.0668 


L6x6-12(F) 

17.3295 

-17.3295 

1.7392 

-1.7392 

-0.4307 

0.4307 

1 3764 

0.1186 


L6x6-12(FR) 

17.4866 

17.4866 

1.7610 

-1.7610 

-0.4357 

0.4357 

1.3905 

0.1160 


L6x6-12(R) 

17.5540 

-17.5540 

1.7647 

-1.7647 

-0.4368 

0.4368 

1.3919 

0.1156 


L6x6-6(R) 

16.7843 

-16.7843 

1.7173 

-1.7173 

-0.4233 

0.4233 

1.3711 

0.1070 

4 

L6x6-6(FR) 

16.9292 

-16.9292 

1.7385 

-1.7385 

-0.4282 

0.4282 

1.3849 

0.1041 


L6x6-6(R) 

16.9970 

-16.9970 

1.7422 

-1.7422 

-0.4292 

0.4202 

1.3863 

0.1038 


Q3x3-12(F) 

17.4836 

-17.4836 

1.7499 

-1.7499 

-0.4342 

0.4342 

1.3844 

0.1225 


Q3x3-12(R) 

17.4963 

-17.4963 

1.7528 

-1.7528 

-0.4343 

0.4343 

1.3846 

0.1229 


Q3x3-6(F) 

16.9386 

-16.9386 

1.7287 

-1.7287 

-0.4268 

0.4268 

1.3789 

0.1109 


Q3x3-6(R) 

16.9462 

-16.9462 

1.7310 

-1.7310 

-0.4268 

0.4268 

1.3791 

0.1113 


EXACT 

1.10 

1.14 

0.119 

0.109 

0.0281 

0.0260 

0.351 

0.0334 


L6x6-12(F) 

70.0028 

-70.0028 

4.0930 

-4.0930 

-1.1682 

1.1682 

4.0684 

0.1266 


L6x6-12(FR) 

71.2805 

-71.2805 

4.1794 

-4.1794 

-1.1919 

1.1910 

4.1496 

0.1036 


L6x6-12(F) 

71.6336 

-71.6336 

4.1940 

-4.1940 

-1.1966 

1.1966 

4.1504 

0.1034 

10 

L6x6-6(F) 

69.5053 

-69.5053 

4.0728 

-4.0728 

-1.1617 

1.1617 

4.0647 

0.1154 


L6x6-6(FR) 

70.7613 

-70.7613 

4.1582 

-4.1582 

-1.1851 

1.1851 

4.1456 

0.0920 


L6x6-6(R) 

71.1144 

-71.1144 

4.1729 

-4.1729 

-1.1898 

1.1898 

4.1164 

0.0920 


Q3x3- 12(F) 

71.4394 

-71.4394 

4.1527 

-4.1527 

-1.1927 

1.1927 

1.1446 

0.1487 


Q3x3-12(R) 

71.5489 

-71.5489 

4.1751 

-4.1751 

-1.1930 

1.1930 

4.1456 

0.1496 


Q3x3-6(F) 

70.9549 

-70.9549 

4.1331 

-4.1331 

-1.1859 

1.1859 

4.1402 

0.1370 


Q3x3-6(R) 

71.0357 

-71.0357 

4.1557 

-4.1557 

-1.1862 

1.1862 

4.1413 

0.1379 


EXACT 

0.725 

0.726 

0.0435 

0.0418 

0.0123 

0.0120 

0.420 

0.0152 


Table 5: Stresses on Cross-Ply. w constant 
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a/h 

Mesh 

CTx 


Txy 

Tst 

Ty* 

20 

L6x6- 12(F) 

0.2421 

-0.2421 

0.1121 

-0.1121 

-0.3444 

0.3444 

0.8048 

0.1276 

L6x6-12(FR) 

0.2552 

-0.2552 

0.1183 

-0.1183 

-0.3633 

0.3633 

0.8406 

0.0066 

L6x6-12(R) 

0.2566 

-0.2566 

0.1189 

-0.1189 

-0.3652 

0.3652 

0.8382 

0.0073 

L6x6-6(F) 

0.2417 

-0.2417 

0.1120 

-0.1120 

-0.3438 

0.3438 

0.8039 

0.1117 

L6x6-6(FR) 

0.2547 

-0.2547 

0.1181 

-0.1181 

-0.3626 

0.3626 

0.8395 

0.0232 

L6x6-6(R) 

0.2561 

-0.2561 

0.1187 

-0.1187 

-0.3645 

0.3645 

0.8394 

0.0239 

Q3x3-12(F) 

0.2558 

-0.2558 

0.1175 

-0.1175 

-0.3637 

0.36-17 

0.8573 

0.2403 

Q3x3-12(R) 

0.2566 

-0.2566 

0.1187 

-0.1187 

-0.3649 

0.3649 

0.8574 

0.2422 

Q3x3-6(F) 

0.2554 

-0.2554 

0.1173 

-0.1173 

-0.3640 

0.3640 

0.8562 

0.2228 

Q3x3-6(R) 

0.2561 

-0.2561 

0.1185 

-0.1185 

-0.3642 

0.3642 

0.8563 

0.2250 

EXACT 

0.2600 

-0.2601 

0.1178 

-0.1195 

-0.3696 

0.3727 

0.8686 

0.2344 

50 

L6x6- 12(F) 

1.1745 

-1.1745 

0.4961 

-0.4961 

-0.1573 

0.1573 

l .5298 


L6x6-12(FR) 

1.5394 

-1.5394 

0.6504 

-0.6504 

-0.2062 

0.2062 

1.8006 


L6x6-12(R) 

1.5482 

-1.5482 

0.6540 

-0.6540 

-0.2074 

0.2074 

1.7984 


L6x6-6(F) 

1.1742 

-1.1742 

0.4960 

-0.4960 

-0.1573 

0.1573 

1.5276 


L6x6-6(FR) 

1.5389 

-1.5389 

0.6502 

-0.6502 

-0.2061 

0.2061 

1.7975 


L6x6-6(R) 

1.5477 

-1.5477 

0.6538 

-0.6538 

-0.2073 

0.2073 

1.7953 


Q3x3- 12(F) 

1.5363 

-1.5363 

0.6432 

-0.6432 

-0.2068 

0.2068 

2.1724 

0.5536 

Q3x3-12(R) 

1.5485 

-1.5485 

0.6534 

-0.6534 

-0.2074 

0.2074 

2.1650 

0.5655 

Q3x3-6(F) 

1.5359 

-1.5359 

0.6429 

-0.6429 

-0.2067 

0.2067 

2.1697 

0.5139 

Q3x3-6(R) 

1.5480 

-1.5480 

0.6533 

-0.6533 

-0.2073 

0.2073 

2.1619 

0.5265 

EXACT 

1.5690 

-1.5691 

0.6458 

-0.6475 

-0.2109 

0.2112 

2.1927 

0.5393 

100 

L6x6- 12(F) 

2.7489 

-2.7489 

1.1443 

-1.1443 

-0.3647 

0.3617 

1.3652 


L6x6-12(FR) 

6.1254 

-6.1254 

2.5499 

-2.5499 

-0.8127 

0.8127 

1.3026 


L6x6-12(R) 

6.1604 

-6.1604 

2.5644 

-2.5644 

-0.8173 

0.8173 

1.2849 

. 

L6x6-6(F) 

2.7489 

-2.7489 

1.1442 

-1.1442 

-0.3647 

0.3647 

1.3625 

. 

L6x6-6(FR) 

6.1249 

-6.1249 

2.5497 

-2.5497 

-0.8126 

0.8126 

1.2966 


L6x6-6(R) 

6.1599 

-6.1599 

2.5642 

-2.5642 

-0.8173 

0.8173 

1.2789 


Q3x3- 12(F) 

6.0577 

-6.0577 

2.5056 

-2.5056 

-0.8107 

0.8107 1 

4.3986 

1.0551 

Q3x3-12(R) 

6.1619 

-6.1619 

2.5642 

-2.5642 

-0.8175 

0.8175 

1.3376 

1.1238 

Q3x3-6(F) 

6.0573 

-6.0573 

2.5052 

-2.5053 

-0.8107 

0.8106 

1.3925 

0.9772 

Q3x3-6(R) 

6.1617 

-6.1614 

2.5639 

-2.5640 

-0.8174 

0.8174 

1.3311 

1.0485 

EXACT 

6.2435 

-6.2436 

2.5308 

-2.5326 

-0.8319 

0.8322 

1.3917 

1.0652 


Table 6: Stresses on Cross-Ply.u) constant 


z 

uTlp 15 

<t> 

1.0000 

Exact 

Variable W 

Constant W 

Exact 

Variable W 

Constant W 

-0.530075 

-0.524338 

-0.537309 

0.000000 

0.000000 

0.000000 

0.9667 

-0.441401 

-0.436096 

-0.454397 

0.009646 

0.009699 

0.002594 

0.9333 

-0.371682 

-0.366761 

-0.390129 

0.018931 

0.019036 

0.004859 

0.9000 

-0.318131 

-0.313540 

-0.341897 

0.027862 

0.028021 

0.006856 

0.9000 

-0.318131 

-0.313540 

-0.341897 

0.027862 

0.028021 

0.006856 

0.7667 

-0.202195 

-0.198911 

-0.240958 

0.023203 

0.023310 

0.006476 

0.6333 

-0.124831 

-0.122843 

-0.163878 

0.019096 

0.19158 

0.006252 

0.6333 

-0.124831 

-0.122843 

-0.163878 

0.019096 

0.19158 

0.006252 

0.5000 

0.036839 

0.037103 

0.000000 

0.015443 

0.015466 

0.006177 

0.3667 

0.199435 

0.198031 

0.163878 

0.012160 

0.012145 

0.006252 

0.3667 

0.199435 

0.198031 

0.163878 

0.012160 

0.012145 

0.006252 

0.2333 

0.273535 

0.270900 

0.240958 

0.009165 

0.009115 

0.006476 

0.1000 

0.361385 

0.357215 

0.341897 

0.006388 

0.006304 

0.006856 

0.1000 

0.361385 

0.357215 

0.341897 

0.006388 

0.006304 

0.006856 

0.0667 

0.405321 

0.400719 

0.390129 

0.004599 

0.004543 

0.004859 

0.0333 

0.465517 

0.460434 

0.454397 

0.002469 

0.002440 

0.002594 

0.0000 

0.544384 

0.538772 

0.537309 

0.000000 

0.000000 

0.000000 


Table 7: Inplane displacement and electrostatic potent i;d 


The laminate is constructed with the upper and . lower layers composed of polyvinylidene 
fluoride, PVDF, oriented at 0 degrees, and with the internal two layers composed of 
a cross-ply of graphite-epoxy oriented at 0, and 90 degrees , respectively [0/90]. The 
geometry of the laminate is shown in figure 18. In this example, two loading cases are 
considered. In the first case, the sinusoidal load is applied to the top surface with the the 
top and bottom surfaces and the vertical edges of the laminate grounded. In the second 
case, te sinusoidal potential is applied to the top surface of the laminate with the bottom 
surface grounded. 

The results for the case l/h = 4 are presented on tables 7 through 9 for the case in 
which the sinusoidal load is applied, and on tables 10 through 12 for the case when the 
sinusoidal potential is applied. 

Semi-Analytic Solutions 

Of primary interest in this section is the basic behaivor of the displacement, stress, 
electrostatic potential, and electric displacement of a laminate with embedded piezoelec- 
tric layers. Both the variable-w and constant-w cases are examined below for the static 
behavior, with the added theory of the independent-w approximation studied for the 
dynamic case. 
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z 

Ox 

T *y 

1.0000 

Exact 

Variable W 

Constant W 

Exact 

Variable W 

Constant W 

10.33504 

10.04129 

10.10107 

-0.659511 

-0.639464 

-0.650112 

0.9667 

8.652445 

8.395070 

8.552297 

-0.567147 

-0.549237 

-0.564983 

0.9333 

7.324507 

7.098908 

7.347798 

-0.489060 

-0.472978 

-0.493142 

0.9000 

6.299222 

6.099411 

6.440052 

-0.423323 

-0.408792 

-0.432799 

0.9000 

0.855399 

0.854236 

0.589416 

-0.372215 

-0.359438 

-0.380547 

0.7667 

0.638399 

0.616228 

0.368611 

-0.193170 

-0.185308 

-0.216892 

0.6333 

0.449534 

0.453350 

0.205646 

-0.067341 

-0.063584 

-0.098115 

0.6333 

1.553872 

1.492791 

1.727943 

-0.067341 

-0.063584 

-0.098115 

0.5000 

-0.213979 

-0.210662 

0.00000 

0.034551 

0.034255 

0.000000 

0.3667 

-1.988776 

-1.919539 

-1.727943 

0.131542 

0.127279 

0.098115 

0.3667 

-0.030955 

-0.037756 

-0.205646 

0.131542 

0.127279 

0.098115 

0.2333 

-0.212196 

-0.193350 

-0.368611 

0.247282 

0.239053 

0.216892 

0.1000 

-0.400353 

-0.394469 

-0.589416 

0.403180 

0.390294 

0.389547 

0.1000 

-6.892361 

-6.686350 

-6.44052 

0.458540 

0.443884 

0.432799 

0.0667 

-7.736448 

-7.505797 

-7.347798 

0.516071 

0.499933 

0.493142 

0.0333 

-8.884439 

-8.623567 

-8.552297 

0.585343 

0.567453 

0.564983 

0.0000 

-10.38138 

-10.08549 

-10.10107 

0.668041 

0.648098 

0.650112 


Table S: Inplane stress distribution 


z 

r v-- 

Dl * 10-“ 

0.9833 

Exact 

Variable W 

Constant W 

Exact 

Variable \V j Const. ml YV 

0.030842 

0.027313 

0.025778 

-0.126731 

-0.120614 

-0.019162 

0.9500 

0.086931 

0.082322 

0.077648 

-0.120156 

-0.114803 

-0.016218 

0.9167 

0.136315 

0.130711 

0.123182 

-0.107478 

-0.102818 

-0.010268 

0.8333 

0.435957 

0.410734 

0.405381 

-0.092683 

-0.092463 

-0.00 7448 

0.7000 

0.750776 

0.726830 

0.736019 

-0.081694 

-0.081487 

-0.004399 

0.5667 

0.831595 

0.817600 

0.827206 

-0.072650 

-0.072465 

-0.001452 

0.4333 

0.849577 

0.835698 

0.827206 

-0.065336 

-0.065183 

0.001452 

0.3000 

0.766156 

0.743298 

0.736019 

-0.059577 

-0.059464 

0.004399 

0.1667 

0.426790 

0.403100 

0.405381 

-0.055236 

-0.055172 

0.00744S 

0.0833 

0.126715 

0.121428 

0.123182 

-0.049454 

-0.053851 

0.010268 

0.0500 

0.080648 

0.076351 

0.077648 

-0.043180 

-0.048213 

0.016218 

0.0167 

0.028558 

0.025335 

0.025778 

-0.039831 

-0.045558 

0.019162 


Table 9: Transverse shear stress and electric displacement 
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z 

r yI 

* IQ' 11 

1.0000 

Exact 

Variable W 

Constant W 

Exact 

Variable VV 

Constant W 

-0.464016 

-0.454333 

-0.274944 

1.000000 

1.000000 

1.000000 

0.9667 

-0.750174 

-0.740068 

-0.776016 

0.985118 

0.985081 

0.984102 

0.9333 

-1.093318 

-1.082760 

-1.229627 

0.971528 

0.971463 

0.969594 

0.9000 

-1.505134 

-1.494176 

-1.661744 

0.959215 

0.959128 

0.956455 

0.9000 

-1.505134 

-1.494176 

-1.661744 

0.959215 

0.959128 

0.956455 

0.7667 

-1.016508 

-1.006818 

-1.105801 

0.770200 

0.769887 

0.767967 

0.6333 

-0.609389 

-0.600110 

-0.660315 

0.599521 

0.599117 

0.597904 

0.6333 

-0.609389 

-0.600110 

-0.660315 

0.599521 

0.599117 

0.597904 

0.5000 

-0.364444 

-0.356267 

-0.399691 

0.443115 

0.442720 

0.442186 

0.3667 

-0.234489 

-0.228002 

-0.268975 

0.297258 

0.296946 

0.297076 

0.3667 

-0.234489 

-0.228002 

-0.268975 

0.297258 

0.296946 

0.297076 

0.2333 

-0.244599 

-0.239503 

-0.298538 

0.158478 

0.158295 

0.159094 

0.1000 

-0.283474 

-0.279194 

-0.363899 

0.023471 

0.023443 

0.024929 

0.1000 

-0.283474 

-0.279194 

-0.363899 

0.023471 

0.023443 

0.024929 

0.0677 

-0.311356 

-0.307364 

-0.371325 

0.015630 

0.015611 

0.016599 

0.0333 

-0.367545 

-0.363802 

-0.351645 

0.007809 

0.007800 

0.008294 

0.0000 

-0.454328 

-0.450824 

-0.304332 

0.000000 

0 . 000000 . 

0.000000 


Table 10: Inplane displacement and electrostatic potential 


Z 

c x * 10 1 

T X y * 10 ^ 

1.0000 

Exact 

Variable VV 

Constant W 

Exact 

Variable VV 

Cons tail L VV 

-0.252476 

-0.238020 

-0.307042 

0.061703 

0.078828 

-0.197721 

0.9667 

-0.168379 

-0.169224 

-0.196349 

-0.508764 

-0.480931 

-0.572105 

0.9333 

-0.074172 

-0.076568 

-0.079646 

-1.072601 

-1.033926 

-1.223532 

0.9000 

0.032380 

0.014740 

-0.017074 

-1.643345 

-1.593534 

-1.783186 

0.9000 

0.017152 

0.018635 

0.026079 

-1.444943 

-1.401145 

-1.567900 

0.7667 

0.012998 

0.012217 

0.017026 

-0.941803 

-0.908786 

-1.007435 

0.6333 

0.009345 

0.007959 

0.010671 

-0.616103 

-0.590772 

-0.656807 

0.6333 

0.067112 

0.065187 

0.072009 

-0.616103 

-0.590772 

-0.656807 

0.5000 

0.041152 

0.039400 

0.044246 

-0.436790 

-0.416920 

-0.469842 

0.3667 

0.027072 

0.025825 

0.030211 

-0.330848 

-0.314971 

-0.363867 

0.3667 

0.005157 

0.004695 

0.005225 

-0.330848 

-0.314971 

-0.363867 

0.2333 

0.004410 

0.004098 

0.005279 

-0.303844 

-0.289910 

-0.346856 

0.1000 

0.003844 

0.004317 

0.006324 

-0.346085 

-0.332338 

-0.410615 

0.1000 

-0.116128 

-0.115099 

-0.117675 

-0.393606 

-0.377970 

-0.466996 

0.0667 

-0.109948 

-0.109214 

-0.115478 

-0.500968 

-0.483674 

-0.562701 

0.0333 

-0.098760 

-0.098185 

-0.118428 

-0.605557 

-0.586605 

-0.595639 

0.0000 

-0.082135 

-0.081793 

-0.127015 

-0.709345 

-0.688722 

-0.566045 


Table 11: Inplane stress distribution 
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z 

Ty, * IQ " 3 

£>, * 10“ l ° 

0.9833 

Exact 

Variable W 

Constant YV 

Exact 

Variable W 

Constant W 

-0.540721 

-0.497105 

-0.568975 

-0.506163 

-0.500145 

-0498431 

0.9500 

-1.510913 

-1.455750 

-1.668638 

-0.462265 

-0.456622 

-0.454952 

0.9167 

-2.336173 

-2.270831 

-2.602503 

-0.418973 

-0.413703 

-0.412115 

0.8333 

-1.757666 

-1.768708 

-2.025278 

-0.376004 

-0.371419 

-0.369941 

0.7000 

-0.398307 

-0.407344 

-0.474639 

-0.339528 

-0.335167 

-0.333779 

0.5667 

0.194760 

0.192967 

0.205508 

-0.311136 

-0.306955 

-0.305625 

0.4333 

0.330157 

0.324822 

0.354593 

-0.290150 

-0.286109 

-0.284803 

0.3000 

0.677038 

0.660715 

0.739205 

-0.276072 

-0.272126 

-0.270814 

0.1667 

1.211837 

1.192111 

1.364158 

-0.268567 

-0.264672 

-0.263323 

0.0833 

1.236155 

1.218844 

1.410369 

-0.266737 

-0.262607. 

-0.261199 

0.0500 

0.733438 

0.724337 

0.841775 

-0.266039 

-0.261917 

-0.260464 

0.0167 

0.242057 

0.241098 

0.279375 

-0.265690 

-0.261570 

-0.260096 


Table 12: Transverse shear stress and electric desplacement 


Quasi-static response 

In this section, the static response of piezoelectric laminates is studied and compared 
with the exact three-dimensional solution. Two different types of loading are considered. 
The first is an applied sinusoidal transverse load of the form 


F(x, y) = fo sin y- sin y- (290) 

where f„ is the peak intensity of the load at the center of the plate. This could simulate 
the sensory characteristics of a laminate by determining the behavior of the electric field 
as a function of loading. The second type of loading simulates the active response of 
a laminate under a sinusoidal surface potential similar in nature to that described in 
Equation 287 except that now f 0 is the peak potential at the plate center. Both of these 
cases have exact solutions [S7] and demonstrate not only the fundamental behavior of the 
laminate but also the accuracy of the discrete- layer approach as a function of the number 
of layers used within the laminate. 

Two geometries are used in this example: a single layer of PZT-4 and a 5-ply hybrid 
laminate. The 5-ply laminate geometry consists of a symmetric [0/90/0] cross-ply of an 
elastic, orthotropic plate with the material properties Cn = 134.9 (all in GPa), C22 = 
14.35, C33 = 14.35, C 12 = 5.156, C 13 = 5.156. C 23 = 7.133, C 44 = 3.606, C 55 = 5.654, Cee 
= 5.654, en/e 0 = 3.5, £22/^0 = £33 jt 0 = 3.0. These three layers are all of equal thickness. 
Two layers of the piezoceramic material PZT-4 of equal thickness are bonded to the upper 
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and lower surfaces of the laminate. The total laminate thickness is defined as h, with the 
thickness of each of the piezoelectric layers taken as O.lh. The plate is square, with the 
length/thickness ratio of the plate specified as 4. 

The boundary conditions along the vertical edges are those consistent with the func- 
tions described in Equation 25. Hence along an edge the transverse and tangential dis- 
placements are zero as are the normal stress and electrostatic potential. The top and 
bottom surfaces of the laminate are grounded for the transverse load case, with the bot- 
tom of the laminate stress-free. For the sinusoidal unit potential, the bottom surface of 
the laminate is grounded with the upper and lower surfaces of the laminate stress-free. 

Three different discretizations are used in the thickness direction of the laminate. Each 
of the laminae are divided into 1, 2, and 4 layers with linear interpolations used within 
each layer. Hence the total laminate is divided into 5, 10, and 20 layers. The in-plane 
functions are selected to coincide exactly with the distributions given in Equation 25, 
which also match the forms for the exact solutions. 

The through-thickness distributions of u, w, (f>, a y , a z ,r yz , and D ; are shown in Figures 
65-80 for the single layer of PZT-4 using a 12 layer discretization. This example clearly 
demonstrates the poor behavior of the constant-w theory and the excellent agreement 
with the variable-w theory. 

The results for the 5-ply laminate in Figures 81-102. The in-plane displacements for 
both applied load an potential contain distinct breaks in slope at the interface locations. 
This is especially true for the applied potential, for which the displacement gradient in 
the PZT layer is very high compared to the other layers. The transverse displacement 
w also has a highly non-uniform behavior over the thickness of the laminate. These 
distributions indicate that even for a relatively simple lamination scheme, the assumption 
of linear global behavior for the displacements or potential would be highly detrimental 
for a plate of this thickness. Also of note are the excellent results obtained using the 
minimal number of discrete-layers through the thickness. Even the 5-laver approximation 
(one layer per laminae) provides excellent results for all field distributions. The stress 
variables in this case are computed at the sub-layer centroids using the constitutive laws 
for the material. The worst agreement for the 5-layer case is the transverse displacement 
for the applied potential loading, for which the values differ from the exact solution by 
about seven percent. This also influences the cr z and r yz stress distributions. All other 
values are in very good agreement for both load cases. 
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The distribution of the normal component of electric displacement is fairly uniform 
through the thickness of the cross-ply by changes dramaticlly within the piezoelectric 
layers because of the electromagnetic coupling. Assuming linear material behavior, this 
behavior could be used to model the laminate response in a sensory fashion for the case 
of the applied loading. 

Dynamic Analysis 

Exact solutions for the free- vibration of simply-supported piezoelectric plates are avail- 
able [89], and provide a good benchmark for comparative purposes. The primary quanti- 
ties of interest are the resonant frequencies and the through-thickness modal distributions 
corresponding to these frequencies. Two sets of boundary conditions are considered in 
this analysis. In the first, the top and bottom surfaces are grounded during the vibra- 
tion. In the second, the electric displacement on these surfaces are zero. These cases are 
termed closed-circuit and open-circuit, respectively. The fundamental in-plane mode with 
m=n=l is the focus of this example. 

A square plate composed of a single layer of the piezoceramic PZT-4 is considered 
first. The length of each side is taken as L x = L y = a. The height is taken as h, and three 
a/h ratios are studied: 4, 10, and 50. The fundamental through- thickness frequency is 
of most interest, and is given as a function of the number of layers used to describe the 
piezoelectric layer. The results are shown in Table 13. The frequencies are represented in 
the tables in terms of the parameter u> h/p. 

A second plate is composed of two dissimilar piezoelectric materials modeled after 
PZT-4 [6S] and PVDF [16]. The densities of the two materials are taken to be the same. 
A three-ply laminate is constructed with the configuration [PZT/PVDF/PZT], with the 
PVDF layer oriented at 0 degrees. The thickness of each PZT layer is 0.25h. Both open 
and closed circuit conditions are considered, with the length/thickness ratios of 4 and 50 
studied. The convergence of the first six modes as a function of a number of sub-layers is 
shown and compared with the exact frequencies in Table 14. For both of these examples, 
it is clear that even a small number of layers yields frequencies accurate well within several 
percent. 

A final example is a 5-ply laminate identical to that considered in the static analysis 
in the previous sub-section. The first six thickness mode frequencies are listed in Table 
15, and are again very accurate with respect to the exact solution. The mode shapes can 
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also be easily calculated, and exhibit trends similar to the static example in terms of the 
convergence. The through-thickness distributions of the displacements and potential for 
the closed-circuit boundary conditions can be easily computed for a/h=4 and a/h=50. 
The plots are indistinguishable from those generated using the exact solution, and are 
given in the next section. 

Exact Solution 

Although both the static and dynamic cases have been investigated, the dynamic case 
is of most interest and can be obtained as a special case of the static solution with u> = 
0. The focus of the results in this report is for the dynamic case. 

In all cases studied, the plates are square with simply supported edges and L x = L y = 
a. The plate surfaces are assumed to be traction-free unless otherwise noted. Although 
other surface conditions can easily be considered along with rectangular plate geometries, 
these were selected for simplicity. Two sets of electric boundary conditions are typically 
applied: homogeneous potential or homogeneous electric displacement. The frequencies 
for these examples are expressed in terms of the frequency parameter *, = wh/^//?, where 
u: is the initial frequency in radians per second. 

The frequencies were determined using bisection with the bounding guesses determined 
two different ways. In the first, the frequency was stepped in increments of one percent of 
the elastic plate frequency neglecting the electromechanical coupling and using classical 
lamination theory [88]. At locations where the determinant changed sign, bisection was 
used to refine the roots to the required accuracy. A second method used as a check was the 
discrete-layer theory developed by the authors [89], which is extremely accurate and gave 
excellent estimates for the exact values of the frequencies. For the most part, these two 
different sets were in very good agreement. The only exceptions existed in the stepping 
method where there could be sign change in determinant and convergence obtained for 
a frequency, but there was no analogous root predicted by the discrete-layer theory. In 
these cases, the resulting exact eigenfunctions did not satisfy the appropriate interface 
and boundary conditions. The source of this anomaly is unknown. Such modes and their 
frequencies are not included in the results that follow. 

The Piezoelectric Single Layer 

A single homogeneous layer of a piezoceramic is considered first. This problem is useful 
to demonstrate the nature of the thickness modes and to partially demonstrate the influ- 
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ence of aspect ratio and the electric boundary conditions. The material properties used 
here, selected to model PZT-4, are shown in the Table A1 in the Appendix along with the 
properties for all materials used in this study. The material used here is material 2. Four 
thickness (a/h) ratios were considered: 1, 4, 10, and 50. The frequencies corresponding 
to the first six thickness modes are shown in Table 16 for both <j >= 0 and D,=0 conditions, 
which are referred to as cases I and II, respectively. It is clear from these results that as 
the aspect ratio decreases (thick plates), the influence of the electric boundary conditions 
becomes more pronounced, with the homogeneous electric displacement providing higher 
frequencies for all cases. 

The mode shapes corresponding to the frequencies of vibration are also of significant 
interest, as these describe the nature of the motion and the extent of electromechanical 
coupling. The three displacement components are normalized with respect to the largest 
value of u, v, or w through the thickness for a particular mode. The potential was also 
normalized with respect to its largest value as well for plotting purposes. A scaling factor 
a was used to denote the relationship between the potential and the displacements. This 
value is the magnitude of the potential at the middle surface of the plate (z=0) divided 
by the u displacement component at the top of the plate (z=h/2). If the potential is zero 
at the mid-surface, the location used for the computed scale factor is at z=0.25h. This is 
noted in the table by placing parentheses around the value. 

The first mode for this laminate is the so-called flexural mode. It is distinguished 
by the symmetric distributions for the transverse displacement and potential and the 
antisymmetric distribution of in-plane displacement. As the thickness ratio increases, 
the in-plane displacements become more linear and the transverse displacement becomes 
more uniform. The second mode is a purely extensional mode with the eigenfunctions 
described by u(z) = -v(z) = 1, w(z) = <f>( z) = 0. There is no electromechanical coupling 
in this mode, and the same frequency is obtained for both case I and case II surface 
conditions. The third mode yields displacement functions that tend toward u(z) = v(z) 
= 1 and w(z) = 0 as the plate becomes thin. The in-plane functions are symmetric and 
the transverse displacement is antisymmetric. The potential distribution changes little 
relative to the aspect ratio. The fourth mode is also purely elastic with no coupling and 
is the first thickness shear mode. This is distinguished by in-plane displacement functions 
given by u(z) = -v(z) = sin 7rz/h and zero transverse displacement and potential. These 
modes demonstrate the deviations in the displacement and electrostatic potential from 
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the kinematic assumptions of some simplified laminate theories (linear u, v, and 4> and 
constant w) as the plate a/h ratio decreases. 

In general, the modes for the two types of electric boundary condition have similar 
features. The features can change depending on the electric boundary conditions. For 
example, for a/h= 1 the second thickness shear mode corresponds to mode 6 (case I) and 
mode 5 (case 2). This mode has features similar to the first thickness-shear mode except 
u(z) = -v(z) = sin 37rz/h. 

Three-ply Symmetric Laminate 

A second example is a laminate composed of two dissimilar materials with a mismatch 
in both elastic and electric properties. The materials, denoted as 1 and 2. are adapted from 
the properties of the transversely isotropic PZT-4 and the orthotropic PVDF. Two aspect 
ratios of 4 and 50 are considered, with both types of electric boundary conditions, I and 
II, considered. Two lamination schemes are studied. The first has the layup of [1/2/2/1], 
and the second [2/1/1/2], where the numbers indicate the piezoelectric material. Each 
layer has equal thickness of 0.25h. The frequency parameters are shown in Table 17. In 
this case, there are no cases for which the elastic and electric fields uncouple. It is also 
not possible to classify these modes as pure shear or extension because of the dissimilar 
materials. 

Plots of the first six through-thickness modes for [1/2/2/1] are shown in Figures 103- 
114 for the two aspect ratios. It is clear from these plots that the influence of the dissimilar 
properties on the displacements and the potential decreases as the aspect ratio increases. 
It is also clear that even for thin laminates the modal potential cannot be accurately 
represented by a simple linear function through the laminate thickness. This requires 
specific attention in constructing approximate solutions to this class ul problem. 

Hybrid Composite Laminate 

In many structural applications, several layers of piezoelectric material are bonded to 
a substructure of elastic composite plies. This type of configuration is considered here, 
with single layers of piezoceramic material PZT (material 2) bonded to the upper and 
lower faces of an elastic, symmetric cross-ply. The 5-ply laminate [PZT/0/90/0/PZT] is 
considered, with the cross-plies formed by a composite material and denoted as material 
3 in Table Al. 

The natural frequencies for case I and II electric boundary conditions are shown in 
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Table 18. The through-thickness displacement and potential functions are presented for 
the first (flexural) mode in Figures 115-116 for a/h=4 and a/h=50. The line definitions 
are given by u (solid), v (dash-dot), w (dash), and the potential (dot). In addition, 
representative plots of the intralaminar (oi, cr 2 , cr 6 ) and interlaminar (cr 3 , o\j, cr 5 ) stress 
components and transverse electric displacement (D 2 ) are also given for the same mode 
in Figures 117-120, respectively, Here the lines are a x (solid), a y (dash), a xy (dot) in 
Figures 117-118 and a z (solid), t xz (dash), T yz (dot) and D 2 (dash-dot) in Figures 119- 
120 for both two thickness ratios. As the plate becomes thin, the in-plane displacements 
tend to the same distribution, with the transverse displacement tending to a constant 
through-thickness value. Also, the in-plane stresses tend to become more linear, with the 
transverse normal stress decreasing in relative magnitude as expected. 

Single-layer PZT Independent-w theory 

The final example of this report contains perhaps the most iinporiaui example oT this 
study, and provides an example of the theory constructed to combat the poor performance 
of the constant-w theory described in earlier sections. By allowing a transverse normal 
strain that is still less computationally intensive than the total variable- w theory, the 
actuation strain can be captured to yield results that are much closer to reality than the 
constant-w case. 

As a simple demonstration, the closed-circuit free-vibration behavior of a single layer 
of PZT-4 is examined using 1) variable-w theory with two layers, 2) u-v-t6 approximations 
using two layers and w-approximation using a single layer, and 3) the exact solution. The 
results of the first 6 frequencies are shown in Table 19. Clearly, the independent-w theory 
yields a higher fundamental mode, but is indicative of the ralatively good accuracy with 
a small number of w-layers. Many of the higher modes (pure extension and those with 
no coupling) exactly capture the frequency of the plate. 

This theory was just completed and debugged code produced at the end of this study. 
The full power and versatility of this theory and element await full exploration. 
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a/h 

4 

10 

50 

N 

Open 

Closed 

Open 

Closed 

Open 

Closed 

1 

110,189 

104,105 

20,485.7 

20,217.4 

848.586 

848.094 

2 

102,599 

100,294 

18,761.5 

18,656.2 

773.651 

773.458 

4 

99,454 

97,915 

18,255.4 

18,180.8 

753.622 

753.525 

8 

98,456 

97,186 

18,122.4 

18,055.7 

748.576 

748.452 

16 

98,311 

96,994 

18,088.9 

18,024.0 

747.299 

747.177 

32 

98,251 

96,946 

18,080.6 

18,016.0 

746.979 

746.859 

64 

98,236 

96,934 

18,078.4 

18,014.0 

746.898 

746.779 

Exact 

98,232 

96,930 

18,077.8 

18,013.4 

746.873 

746.752 


Table 13 Convergence of fundamental frequency for single-layer of PZT-4. 



Mode 

a/h=4: C 

' 1 

2 

3 

4 

5 

6 

N=4 

74.7732 

197.353 

329.301 

342.387 

447.653 

549.479 

N=8 

72.8253 

195.402 

311.843 

338.628 

432.976 

534.030 

N=16 

72.3376 

194.920 

307.606 

337.504 

426.810 

530.343 

N=32 

72.2152 

194.800 

306. 55S 

337.207 

425.160 

529.432 

N=64 

72.1846 

194.770 

306.296 

337.132 

424.742 

529.205 

Exact 

72.1744 

194.760 

306.209 

337.107 

424.602 

529.129 

a/h=4: 0 

1 

2 

3 

4 

5 

6 

N=4 

74.8023 

197.483 

329.691 

342.481 

447.733 

550.079 

N=8 

72.8451 

195.525 

312.190 

338.718 

433.049 

534.493 

N=16 

72.3554 

195.042 

307.940 

337.593 

426.875 

530.769 

N=32 

72.2325 

194.921 

306.SSS 

337.296 

425.223 

529.849 

N=64 

72.2017 

194.891 

306.626 

337.221 

424.801 

529.620 

Exact 

72.1915 

194.881 

306.539 

337.196 

424.664 

529.543 

a/h=50: C 

1 

2 

3 

4 

5 

6 

N=4 

0.637707 

16.432S 

28.5365 

292.035 

379.322 

400.300 

N=8 

0.634494 

16.4315 

28.5356 

273.945 

359.391 

379.362 

N=16 

0.633687 

16.4312 

28.5353 

269.563 

354.642 

371.982 

N=32 

0.633485 

16.4311 

28.5353 

268.479 

353.469 

370.048 

N=64 

0.633434 

16.4311 

2S.5353 

268.208 

353.177 

369.559 

Exact 

0.633417 

16.4311 

28.5352 

268.118 

353.079 

369.396 

a/h=50: 0 

1 

2 

3 

4 

5 

6 

N=4 

0.637786 

16.4426 

28.5566 

295.866 

390.811 

400.301 

N=8 

0.634566 

16.4413 

28.5556 

277.248 

369.149 

379.363 

N=16 

0.633757 

16.4410 

28.5554 

272.717 

363.958 

371.982 

N=32 

0.633555 

16.4409 

28.5554 

271.595 

362.674 

370.04S 

N=64 

0.633506 

16.4409 

2S.5553 

271.315 

362.354 

369.560 

Exact 

0.6334S7 

16.4409 

28.5553 

271.222 

362.248 

369.396 


Table 14 Convergence of frequencies for 3-ply piezoelectric laminate. 
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Mode 

a/h=4: C 

i 

2 

3 

4 

5 

6 

N=5 

57.2531 

194.840 

255.648 

282.168 

368.461 

389.525 

N=10 

57.1249 

192.190 

252.024 

276.853 

364.218 

383.352 

N=20 

57.0875 

191.524 

251.085 

275.425 

362.939 

381.629 

N=32 

57.0855 

191.350 

250.845 

275.040 

362.592 

381.171 

N=64 

57.0773 

191.313 

250.786 

274.969 

362.516 

381.072 

Exact 

57.0745 

191.301 

250.769 

274.941 

362.522 

381.049 

a/h=4: 0 

1 

2 

3 

4 

5 

6 

N=5 

57.2707 

194.843 

255.648 

2S2.16S 

36S.505 

389.534 

N=10 

57.1403 

192.192 

252.025 

276.853 

364.252 

383.364 

N=20 

57.1023 

191.526 

251.086 

275.425 

362.971 

381.641 

N=32 

57.1005 

191.353 

250.845 

275.040 

362.623 

381.184 

N=64 

57.0921 

191.316 

250.7S6 

274.969 

362.548 

381.084 

Exact 

57.0S93 

191.304 

250.770 

274.941 

362.522 

381.049 

a/h=50: C 

1 

2 

3 

4 

5 

6 

N=5 

0.619025 

15.6835 

21.4947 

212.S11 

214.690 

334.953 

N=10 

0.61834S 

15.6820 

21.4933 

210.561 

211.596 

379.943 

N=20 

0.61S175 

15.6817 

21.4929 

209.925 

210.791 

378.575 

N=32 

0.61S156 

15.6S16 

21.492S 

209.754 

210.582 

378.207 

N=64 

0.61S127 

15.6816 

21.492S 

209.71S 

210.538 

378.132 

Exact 

0.61S11S 

15.6S16 

21.4928 

209.704 

210.522 

378.104 

a/h=50: 0 

1 

2 

3 

4 

5 

6 

N=5 

0.61903S 

15.6S35 

21.4949 

212.827 

214.736 

384.953 

N=10 

0.618351 

15.6821 

21.4935 

210.56S 

211.645 

379.944 

N=20 

0.618179 

15.6817 

21.4931 

209.929 

210.841 

378.575 

N=32 

0.618160 

15.6816 

21.4930 

209.758 

210.632 

373.207 

N=64 

0.618141 

15.6816 

21.4930 

209.721 

210.589 

378.133 

Exact 

0.61S120 

15.6816 

21.4930 

209.707 

210.573 

378.105 


Table 15 Convergence of frequencies for 5-ply hybrid laminate. 
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I 


II 

a/h=l 

7 

a 

7 

a 

i 

713,061 

-1.744e+9 

724,602 

-1.675e+9 

2 

777,021 

0 

777,021 

0 

3 

889,902 

(-1.515e+ll) 

912,912 

(-1.49Se+9) 

4 

925,431 

0 

925,431 

0 

5 

1,243,819 

-1.534e+9 

1,270.594 

0 

6 

1,270,594 

0 

1.293,504 

(-2.57e+10) 

a/h=4 

7 

a 

7 

a 

1 

96,929.9 

-4.87e+8 

98,231.7 

-7.78e+S 

2 

194,255 

0 

194,255 

0 

3 

327,663 

(2.276e+7) 

355,110 

(4.57e+9) 

4 

538,885 

0 

538,885 

0 

5 

609,186 

- 1 .OSe-pS 

690,767 

2.55e+9 

6 

958,922 

(-2.66e+S) 

960,103 

(-3.4Se+S) 

a/h=10 

7 

a 

7 

a 

1 

18,013.4 

-1.78e+8 

IS, 077. S 

-3.40e+8 

2 

77,702.1 

0 

77,702.1 

0 

3 

133,695 

(L27e+6) 

145,221 

(1.74e+S) 

4 

508,625 

0 

508,625 

0 

5 

522,320 

-4.87e-r7 

604,752 

6.26e+9 

6 

988,021 

(-1.76e-i-S) 

990,953 

(-2.8Se+S) 

a/h=50 

7 

a 

7 

a 

1 

746.752 

-3.48e + 1 

746.873 

-7.78e+S 

2 

15,540.4 

0 

15,540.4 

0 

3 

26,828.0 

(9.94e+3) 

29,153.3 

(3.48e+ 1 ) 

4 

502,895 

0 

502,895 

0 

5 

503,469 

-1.01e+7 

586,240 

3.10(:*+ 10 

6 

1,004,344 

(-4.56e+7) 

1,004,612 

(-7.21e+7) 


Table 16. Frequency parameters for single-ply piezoelectric layer for differing electric 

surface conditions. 


SS 



Frequency parameter 7 


a/h=4 a/h=50 


Mode I II I II 


1 

72,174.4 

72,191.5 

633.417 

633.487 

2 

194,760 

194,881 

16,431.1 

16,440.9 

3 

216,505 

216,505 

17,320.4 

17,320.4 

4 

306,209 

306,539 

28,535.2 

28,555.3 

5 

337,107 

337,196 

268,118 

271,222 

6 

424,602 

424,664 

353,079 

362,248 


17a. [ 1 / 2 / 2 / 1 ] lamination scheme. 



Frequency parameter 7 


a/h=4 

a/h 

=50 

Mode 

I 

II 

I 

II 

1 

58,248.7 

58,354.0 

725.219 

725.241 

2 

192,408 

192,436 

16,430.2 

16,438.8 

3 

271,757 

271,758 

28,535.7 

28,555.1 

4 

329,584 

329,593 

159,732 

159,865 

5 

363,048 

364,072 

226,218 

226,643 

6 

406,665 

407,771 

353, 3S6 

363,810 


17b. [2/ 1/1/2] lamination scheme. 


Table 17. Frequency parameters for laminate with dissimilar piezoelectric layers. 
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Frequency parameter 7 


a/h=4 

a/h 

=50 

Mode 

I 

II 

1 

II 

1 

57,074.5 

57,089.3 

618.118 

618.120 

2 

191,301 

191,304 

15,681.6 

15,681.6 

3 

250,769 

250,770 

21,492.8 

21,493.0 

4 

274,941 

274,941 

209,704 

209,707 

5 

362,492 

362,522 

210,522 

210,573 

6 

381,036 

381,049 

378,104 

378,105 


Table 18. Frequency parameters for 5-ply elastic/piezoelectric laminate. 


Frequency parameter 7 
Mode Variable w Linear w Exact 


1 

100.293.4 

104145.2 

96,929.9 

2 

194,255 

194,255 

194,255 

3 

328,560 

328,560 

327,663 

4 

587.312 

587,312 

508,625 

5 

673.851 

706,727 

522,320 

6 

1,125,404 

1,125,404 

9SS,021 


Table 19. Comparison of frequencies for variable w, linear w, and exact, single layer 

PZT-4. 




90 


8. SUMMARY AND FUTURE WORK 


Developed Algorithms 

The following algorithms have been developed as part of this research to aid in the 
simulation of active tip clearance control using adaptive composites: 

• Semi-analytic solutions using variable, constant, and independent approximations 
for out-of-plane transverse displacements for laminated piezoelectric plates. 

• Finite element approximations using constant and variable variations. Implementa- 
tion of independent-w theory is straightforward. Both static and dynamic problems 
can be modeled. 

• Exact solutions of simply-supported laminated piezoelectric plates: static and dy- 
namic behavior. 

• Discrete-layer shell elements for arbitrary geometry. 

Deliverables to NASA-Lewis are the source code for the plate and becim elements, and 
have been delivered to NASA via internet. All other source codes for the semi-analytic 
and exact solutions is also available. 

Summary of Results 

Inspection of the results included here and as part of this study, the primary thrusts 
of the working algorithms are: 

• Excellent accuracy for variable-w theory and promising results for the independent- 
w theory. 

• As currently formulated, the constant-w theory is not suitable for coupled problems 
in electroelasticity because of the inability to model transverse normal strain. For 
elastic problems, however, this methodology yields very acceptable results. 

• Exact solutions developed will provide extremely useful benchmarks for the enclosed 
theories and those developed by others. 


91 


• Flat-plate results for plate element are good and encouraging. The application to 
an assemlage of flat plates to approximate a curved surface has not been made, but 
can be with the inclusion of developed transformation matrices. 

• Elements developed here can be incorporated into general purpose finite element 
code to simulate the control process for a wide array of problems. 

• Results of shell element are promising, but this application has not been extensively 
tested because of lack of comparative solutions available in literature and elsewhere. 

Future Work 

The potential applications of the elements developed here are numerous. Before this 
occurs, additional testing of both the plate and shell elements would be prudent and 
is necessary for the shell element. What follows are suggestions for future work by the 
present investigators or others: 

• Complete implementation and testing of plate element with independent-w theory. 

• Use of rotated plate elements as collection of flat plates to model curved shell. 

• Development of exact solution for laminated piezoelectric cylindrical shell for pur- 
poses of benchmarking general shell element. 

• Further testing of shell element and comparison of each discrete-layer theory. 

• Application to active tip clearance control and other problems requiring adaptive 
composite laminates. 


I 
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Figure 1: Coordinate axes for typical piezoelectric layer. 



Figure 2: Laminated piezoelectric plate and possible deformation. 
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Non-layer- wise description for either u, v, w or <f>. 



Figure 3: Description of pseudo- layer within laminate gem in •try. 



Figure 4: Curvilinear coordinates in a general shell. 
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Layer-wise description for u, v and <f>. 




Figure 5: Through-thickness approximations for discrete- layer theories. 


layer-wise description for u , v, and <j>. 



Figure 6: Assumed form for independent-\v approximation. 


[Kn] = 



a=l 

b=l* 


a^=m 

b=l 



Figure 7: Structure of sub-matrices. 
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There are 32 nodes required for 
this 3-layer shell. 


Figure S: Discrete-layer element for general shells. 



Figure 9: Continuity requirements at through-thickness locations. 
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Figure 10: Shell element with different layers/coordinate systems. 
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Figure 11: Form of approximations for general shell element. 
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Uj(N+l), V j(N+l), W j(N+l), <t> j(N+l) 
UjN , V jN , WjN , <I> jN 


Uj3 ,V jj ,Wb , 4> fl 
Uj2 ,V j2 ,Wj2 . <1> J2 
Uj! ,Vj! ,Wji , <D j! 



Figure 12: Geometry and layer numbering. 

W varying linearly through thickness 


Uj(N+l), V j(N+l), (t> j(N+l) 
UjN , V jN , <t> jN 


Uj3 ,Vj3 , <D j3 

Uj2 , V j2 , 4> j2 
Uji , V j! , <D j! 



Figure 13: Geometry and layer numbering. W constant through thickness 


x, y ,z : global system of coordinates 



x\ y\ z* : local system of coordinates 
{} : angle of rotation between x and x’ 


Figure 14: local and global coordinates. 


q=q o sin(xcx/a) sin(icy/b) 




0 

m 

90 

0 


t- 4 H 



Figure 15: Geometry of Rectangular Plate 
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y 



Figure 16: Borydary Conditions for the Quarter Plate 


q = (jj sin (n x/a)sin (jt y/b) 


or 


= <t» 0 sin (7t x/a)sin (n y/b) 


/fTTTK 




Figure 17: Geometry of PVDF Square Plate 
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q = q 0 sin( 7Z x/1 )sin( n y/1 ) 




Figure 18: Geometry of PVDF/GRAPHITE/EPOXY/PVDF Square Plate 
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Q0C7 

.... LEk6" 1 2 _ r, variable v 
16*6-1 2" cone tan l w 



Figure 19: Displacement in x-direction 


DtflCT 

.... L6x6-12“f\ variable w 
L6*6~12-r, cone tan l w 



-11.0 >1X4 -U HJ -1.0 0.0 U 0.0 0.0 10.0 IS.O^ 

Oweplo c eeenC wn ydlr-ecUon, V e|Q 


Figure 20: Displacement in y-direction 
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— doct 

.... L5*6-12-r, voriabls w 
L6 k 6“ 12 _ r, constant w 



Figure 21: Displacement in z-directiou 


— CTO 

.... L6*6”12-r, voriobls w 
L6n6 - 12 _ r, constant w 



Figure 22: Normal Stress Distribution on x-directiou 
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CXflCT 

.... L6*6-12(n # vorvobl* w 
LSxfi - 1 2 ( r 1 # coritLort « 



Figure 23: Normal Stress Distribution on y-directiuii 


— OtfCT 

.... LSxfi - l2 - r,vorv obi • m 
L6xfi“l2-r , constant » 



Figure 24: In-Plane Stress Distribution 
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DCflCT 

.... Lfix6-12-r, varLabl* w 
L6*6~ 12 _ r, constant w 



Figure 25: Out-plane Stress Distribution 


DCflCT 

.... LS«6-12"F, vorlabls w 
L6*6-12 _ r, constant « 



Figure 26: Out-plane Stress Distribution 


CXACT 

.... L6x6-12~r, voriabla w 
LB*6 _ 12~r # conatanC w 



Figure 27: Displacement in x-direction, applied load ,a/h = 1U. PVDF 


— exfcr 

.... Lfix6*12 _ r # vorvobL* w 
L 6 * 6 “ 12 -r, conatonl w 



Figure 28: Displacement in y-direction with applied load, a/h = 10, PVDF 
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Figure 29: Displacement in z-direction, applied load,a//i = 1U.PVDF 
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Figure 30: Potential Distribution, applied load, a/h = 10.PVDF 
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Figure 31: Normal stress distribution on x-direction, applied load.w/// = 10,PVDF 
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Figure 32: Normal stress distribution on y-direction, applied load ,afh = 10.PVDF 
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Figure 33: Normal stress distribution on z-direction, applied load .a/li = 10,PVDF 



Figure 34: In-plane stress distribution, applied load,«//? = I0.PVDF 
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Figure 35: Out-plane stress distribution, applied load,a//i = 10,PVDF 
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Figure 36: Out-plane stress distribution, applied load, a/li = 1U.PVDF 



Figure 37: Electric displacement distribution, applied load, a/li = 10,PVDF 
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Figure 3S: Displacement in x-direction, applied potential, a/h = LU.PVDF 
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Figure 39: Displacement in y-direction, applied potential, a/h = 10,PVDF 



Figure 40: Displacement in z-direction, applied potential, a/h = 10,PVDF 
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Figure 41: Potential distribution, applied potential, a/h = I0.PVDF 
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Figure 42: Normal stress distribution on x- direction .applied potential 

a/h = 10,PVDF 
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DCflCT 
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Figure 43: Normal stress distribution on y-direttion, applied potential 

a/h = 10,PVDF 



Figure 44: Normal stress distribution on z-direction, applied potential 

a/h = 10,PVDF 
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Figure 45: In-plane stress distribution. applied potential, a/lr = 10,PVDF 



Figure 46: Out-plane stress distribution, applied potential, a/h = 10,PVDF 
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Figure 47: Out-plane stress distribution, applied potential, u/li = 10,PVDF 



Figure 48: Electric displacement distribution, applied potential, a/h = 10,PVDF 
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Figure 49: Displacement in x-direction, applied load,a//i = 10.PZT4 



Figure 50: Displacement in z-direction, applied load,a//i = 10.PZT4 
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Figure 51: Potential distribution, applied load, a/li = 1U.PXT4 
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Figure 52: Normal stress distribution on x-direction, applied load ,u//x — 10,PZ14 
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Figure 53: Normal stress distribution on z-direction, applied \o<u\.<i/h = 
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Figure 54: In-plane stress distribution, applied load ,a/h = 1U.PZT4 
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Figure 55: Out-plane stress distribution, applied load, a//i = 10,PZT4 
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Figure 56: Electric displacement distribution, applied load .<//// = 10.PZT4 
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Figure 57: Displacement in x-direction, applied potential, a/li = 10,PZT4 
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Figure 58: Displacement in z-direction, applied potential, a/h = 10,PZT4 
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Figure 59: Potential distribution, applied potential, a/h = 10.PZT4 
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Figure 60: Normal stress distribution on x-direction, applied potential 

a/h = 10,PZT4 
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Figure 61: Normal stress distribution on z-direction, applied potential 

a/h = 10,PZT4 



Figure 62: In-plane stress distribution, applied potential, a/h = 10,PZT4 
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Figure 63: Out-plane stress distribution, applied potential, a/U = 10.PZT4 



Figure 64: Electric displacement distribution, applied potential, a/ li = 10,PZT4 
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Figure 65: u-displacement for single layer PZT-4 (load). 



Figure 66: w-displacement for single layer PZT-4 (load). 
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Figure 67: Potential for single layer PZT-4 (load). 



Figure 68: In-plane normal stress for single layer PZT-4 (load). 
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ure 69: Out-of-plane normal stress for single layer PZT-d (load). 
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Figure 70: Out-of-plane shear stress for single layer PZT-4 (load). 
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Figure 71: In-plane shear stress for single layer PZT--1 (load). 



Figure 72: Normal electric displacement for single layer PZT--I (load). 
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Figure 73: u- displacement for single layer PZT-4 (potential). 



Figure 74: w-displacement for single layer PZT-4 (potential). 
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Figure 75: Potential for single layer PZT-4 (potential). 



Figure 76: In-plane normal stress tor single layer PZ1-4 (potential). 
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Figure 77: Out-of-plane normal stress for single lai'er PZT-I (potential). 



Figure 7S: Out-of-plane shear stress for single layer FZT-l (potential). 


133 



Figure 79: In-plane shear stress for single layer PZT-4 (potential). 



Figure SO: Normal electric displacement for single layer PZT-I (potential). 
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Figure 81: u-displacement for 5-ply, applied load. 
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Figure 82: v-displacement for 5-ply, applied load. 
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Figure 83: w-displacement for 5-ply, applied load. 



Figure 84: Potential for 5-ply, applied load. 
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Figure 85: Normal (x) stress, 5-ply, applied load. 



Figure 86: Normal (y) stress, 5-ply, applied load. 


137 



Figure S7: Normal (z) stress, 5-plv, applied load. 



Figure SS: Shear (xz) stress, 5-ply, applied load. 
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Figure 89: Shear (yz) stress, 5-ply, applied load. 



Figure 90: Shear (xy) stress, 5-ply, applied load. 
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Figure 91: Normal electric displacement, 5-ply. applied load. 



Figure 92: u-displacement for 5-ply, applied potential. 
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Figure 93: v-displacement for 5-ply, applied potent ial. 
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Figure 94: 


w-displacement for 5-ply, applied potential. 
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Figure 95: Potential for 5-ply, applied potential. 



Figure 96: Normal (x) stress, 5-ply, applied potential. 
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Figure 97: Normal (y) stress, 5-ply, applied potential. 



Figure 9S: Normal (z) stress, 5-ply, applied potential. 
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Figure 99: Shear (xz) stress, 5-ply, applied potential. 



Figure 100: Shear (yz) stress, 5-ply, applied potential. 
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Figure 101: Shear (xy) stress, 5-ply, applied potential. 



Figure 102: Normal electric displacement, 5-ply, applied potential. 
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Figure 103: Mode 1, 3-ply, a/h=4. 
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Figure 104: Mode 1, 3-ply, a/h=50. 
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Figure 105: Mode 2, 3-ply, a/h=4. 



Figure 106: Mode 2, 3-ply, a/h=50. 
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Figure 107: Mode 3, 3-ply, a/h=4. 





Figure 108: Mode 3, 3-ply, a/h=50. 
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Figure 109: Mode 4, 3-ply, a/h=4. 



Figure 110: Mode 4, 3-ply, a/h=50. 
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Figure 111: Mode 5, 3-ply, a/h=4. 
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Figure 112: Mode 5, 3-ply, a/h=50. 
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Figure 113: Mode 6, 3-ply, a/h=4. 



Figure 114: Mode 6, 3-ply, a/h=50. 
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Figure 115: Displacement and potential, mode 1, 5-ply, ;i/li=-|. 



Figure 116: Displacement and potential, mode 1, 5-ply, a/!i=50. 
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Figure 117: In-plane stresses, mode 1, 5-ply, a/li= 1. 



Figure 118: In-plane stresses, mode 1, 5-ply, a/h=50. 
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Figure 119: Out-of-plane stresses, mode 1, 5-ply, a/li= 1. 



Figure 120: Out-of-plane stresses, mode 1, 5-ply, a/li =•'>(). 
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APPENDIX 




Property 12 3 


E x (GPa) 

237.0 

81.3 

132.38 

E 2 

23.2 

81.3 

10.756 

£3 

10.5 

64.5 

10.756 

1/12 

0.154 

0.329 

0.24 

^13 

0.178 

0.432 

0.24 

^23 

0.177 

0.432 

0.49 

G. 44 

2.15 

25.6 

3.606 

G 55 

4.4 

25.6 

5.6537 

G 66 

6.43 

30.6 

5.6537 

e 24 (C/m 2 ) 

-0.01 

12.72 

0 

^31 

-0.13 

-5.20 

0 

632 

-0.14 

-5.20 

0 

633 

-0.2S 

15.08 

0 

ILL 

co 

12.5 

1475 

3.5 

122 

cn 

11.98 

1475 

3.0 

I22 

in 

11.98 

1300 

3.0 


Table Al. Elastic, piezoelectric, and dielectric properties of piezoelectric materials. 
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Definitions of coeeficients used in exact solution: 



(291) 


B — An C44 + An C44 C33 €33 — 2 A13 C44 A14 633 + C55 A24C33 — C55 C44 C33 A44 -{ 292 ) 


C — —2 A12 A23 ^13 ^33 — A u M2 ^33 — C55 A23A44 — C55 A24 A33 + C55 C44 A33 A44 -f{ 293 ) 
^ 12 e 33 “ 2 C55 A22 A34 e33 — C55 C44 A34 + 2 C55 A23 A24 A34 — An A22 C33 633 — 


A12C33 633 + 2 A12 A24 A14 C33 — 2 A 12 A24 A13 633 — An A24C33 — 2 An A 23 Ah 633 + 

2 An A23 A24 633 + AJ3A24 + 2 A13 C44 Ah A34 — 2 An C44 A34 633 + 2 A13 A22 A H e33 + 

^ 14^23 + AJ3A22 ^33 + An A23C33 — 2 A13 A24 A14 A23 “ AJ4A22 ^33 “ A14C44 A33 

D = —An C44 A33 A44 — An A22 C33 A44 + 2 An A22 A34 633 + An A22 A33 £33 + C 55 A22 A34 -{29*1 ) 
C55 A22 A33 A44 + A13A22 A44 — 2 A12 A24 A14 A33 + 2 A12 A24 A13 A34 + 2 A12 A23 Ah A34 — 

2 A12 A23 ^13 A44 - 2 A^ 2 A34 e 33 — A12A33 C33 + A 2 i2 C-&t A 44 + A n A 24 A33 - 
2 An A23 A24 A34 + An -^ 23^44 + An C44 A34 + ^14^22 ^33 — 2 A13 A22 Au A34 




E — An A22 -^33 A44 + A^A^ — An A22 A\ 4 — ^12^33 A44 


( 295 ) 



( 296 ) 


D2 — — C44 Az3 £33 + C44 C33 A44 — 2 C44 A34 633 + 2 A 23 A24 £33 — ( 297 ) 

A22 633 + ^ 23 e 33 — ^22 C33 £33 — ^ j 


165 


D$ = A22 ^33 ^33 — ^22 C33 ^44 — 2 A23 A24 A34 + 2 A22 A34 ^33 + ( 29 S) 

^ 23-^44 “ C44 A33 A44 + ^ 24^33 + C44 ^34 

D4 = A22 A33 A44 ^22 ^34 ( 299 ) 

fll = A\2 (C33 €33 + — (A23 C33 + A24 C33) A13 + ( — A23 ^33 + ^24 ^33) ^14 ( 300 ) 

f\2 = — A\2 (A33 633 — C33 A44 + 2 A34 633) — (A23 A 44 — A24 A34) A\3 + (A23 A34 — A24 ^33) /li( r 5 ()l ) 

/13 = —A12 (A33 A44 — ^4.34^ ( 302 ) 

/21 = “C44 £33 A\3 — C44 633 A14 ( 303 ) 

/22 — (^22 ^33 C44 A44 + A24) A\3 — A\2 ( A23 ^33 + ^24 ^ 33 ) + ( 304 ) 

(-^22 £33 + C44 ^34 “ A23 ^24) A\4 

/23 = — A\2 ( A23 A44 — A24 A34) — A22 A34 A\4 + A22 A44 A 13 ( 305 ) 

/31 = — ^13 ^44 ^33 + ^14 C44 <^33 ( 306 ) 

$32 = A\3 (-A22 ^33 “1" C44 A34 — A23 A24) + A\2 (“^23 e 33 + ^24 C33) + (307) 

A\4 (”'^22 C33 “ C44 ^33 + ^23 ) 

/33 = A\2 (A23 A 34 — A24 ^33) 4 " A\4 A22 A33 A13 A22 A34 ( 30 b) 
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